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ABSTRACT 



Aims. We present the SxAGGER-grid, a comprehensive grid of time-dependent, three-dimensional (3D), hydrodynamic model atmo- 
spheres for late-type stars with realistic treatment of radiative transfer, covering a wide range in stellar parameters. This grid of 3D 
models is intended for various applications besides studies of stellar convection and atmospheres per se, including stellar parameter 
determination, stellar spectroscopy and abundance analysis, asteroseismology, calibration of stellar evolution models, interferometry, 
and extrasolar planet search. In this introductory paper, we describe the methods we applied for the computation of the grid and 
discuss the general properties of the 3D models as well as of their temporal and spatial averages (here denoted (3D> models). 
Methods. All our models were generated with the SxAGGER-code, using realistic input physics for the equation of state (EOS) and for 
continuous and line opacities. Our ~ 220 grid models range in effective temperature, T^ff, from 4000 to 7000K in steps of 500 K, in 
surface gravity, logg, from 1.5 to 5.0 in steps of 0.5 dex, and metallicity, [Fe/H], from -4.0 to +0.5 in steps of 0.5 and l.Odex. 
Results. We find a tight scaling relation between the vertical velocity and the surface entropy jump, which itself correlates with the 
constant entropy value of the adiabatic convection zone. The range in intensity contrast is enhanced at lower metallicity. The granule 
size correlates closely with the pressure scale height sampled at the depth of maximum velocity. We compare the (3D) models with 
currently widely applied one-dimensional (ID) atmosphere models, as well as with theoretical ID hydrostatic models generated with 
the same EOS and opacity tables as the 3D models, in order to isolate the effects of using self-consistent and hydrodynamic modeling 
of convection, rather than the classical mixing length theory (MLT) approach. For the first time, we are able to quantify systematically 
over a broad range of stellar parameters the uncertainties of ID models arising from the simplified treatment of physics, in particular 
convective energy transport. In agreement with previous findings, we find that the differences can be rather significant, especially for 
metal-poor stars. 

Key words. convection - hydrodynamics - radiative transfer - stars: abundances - stars: atmospheres - stars: fundamental 
parameters - stars: general- stars: late-type - stars: solar-type 

1. Introduction for the interaction between radiative and convective energy 

. . ' transport at the optical surface . 

^ The primary source of information for stellar objects is the light 

• they emit, which carries information about the physical condi- 

/\ tions at its origin. However, in order to interpret the informa- The first realistic grids of line-blanketed atmosphere 

H tion correctly, one first needs either theoretical or semi-empirical models for late-type stars appeared with the publication of 

-- -models of the atmospheric layers at the surface of stars from MARCS (Gustafsson et al. 1975, 2008) and ATLAS models 

where the stellar radiation escapes. Therefore, models of stellar (Kurucz 1979; Castelli & Kurucz 2004). Subsequendy, other 

atmospheres are essential for much of contemporary asti-onomy. one-dimensional (ID) atmosphere codes, e.g. PHOENIX 

In the case of late-type stars, the theoretical modeling of (Hauschildt et al. 1999) and MAFAGS (Grupp 2004), were 

stellai- atmospheres is complicated by the presence of convective developed to model the atmospheres of stars. In general, 

motions and turbulent flows as well as of magnetic fields in their these theoretical ID atmosphere models assume hydrostatic 

envelopes (see review by Nordlund et al. 2009, and references equilibrium, flux constancy, and local thermodynamic equilib- 

therein). In particular, convection can significantly affect both num (LTE). For the modeling of convective energy transport, 

the atmospheric stratification and emergent spectral energy they commonly employ the mixing-length theory (MLT, see 

distribution in these stars. Hence, in order to correctly represent Bohm-Vitense 1958), which is characterized by several free 

the temperature stratifications in the outer layers of stai's, from parameters, the most commonly known being the mixing-length 

where the stellai- light escapes, it is vital to accurately account L, or equivalently, the parameter q-mlt = L/Hp. Alternatively, 

some relatives thereof are available, such as the full turbulence 

Send offprint requests to: magic@mpa-garching.mpg.de spectrum (FTS) theory by Canuto & Mazzitelli (1991), which 
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itself also has a free parameter. The values of these free parame- 
ters are not known from first principles and need to be calibrated 
based on observations or simulations. The mixing-length theory 
has in total four free parameters (see Bohm-Vitense 1958; 
Henyey et al. 1965; Mihalas 1970). These free parameters 
can be calibrated based on their effect on synthetic spectra, 
but usually only qtmlt is calibrated based on the reproduction 
of selected lines (Fuhrmann et al. 1993; Barklem et al. 2002; 
Smalley & Kupka 2003). Moreover, the free mixing length is 
calibrated in stellar evolutionary calculations by matching the 
observed luminosity and radius of the Sun at its current age 
(e.g. Magic et al. 2010). To construct simple yet reaUstic ID 
models of convection is rather difficult, in particular convective 
overshooting beyond the classical Schwarzschild instability cri- 
terion is normally not considered in ID atmospheric modeling. 
Attempts have been made at including its effects in ID model 
atmospheres albeit with only limited success (CastelU et al. 
1997). 

The first numerical ID model stellar atmosphere codes 
usually assumed a plane-parallel geometry for the atmospheric 
stratification. This was later improved upon by changing to 
a spherical symmetry, leading to lower temperatures in the 
upper layers, in particular for giant stars, due to the dilution 
of the radiation field with increasing radial distance, which 
can cover a significant fraction of the stellar radius at low 
logg (see Gustafsson et al. 2008). Initially line blanketing was 
included by means of opacity distribution functions (ODFs, 
Gustafsson et al. 1975) with a few hundred ODFs covering the 
entire spectrum, eventually replaced by opacity sampling (OS) 
including thousands of wavelength points (Johnson & Krupp 
1976). Nowadays, thousands of ODFs or hundreds of thousands 
of OS wavelengths are used. Despite such high resolution in 
wavelength, the computational costs for ID atmosphere models 
are currently quite small, at least for LTE models. Large, ho- 
mogeneous grids of atmospheres with up to ~ 10^ models exist 
(Gustafsson et al. 2008; Cassisi et al. 2004; Hauschildt et al. 
1999), covering a wide range of stellar atmosphere parameters 
(Teff, log^, and [Fe/H]). 

Even though the ID atmosphere models are based on nu- 
merous simplifications, they have demonstrated high predictive 
capabihties owing to major improvements in the atomic and 
molecular data (e.g. line lists by Kurucz (1993) or VALD by 
Piskunov et al. (1995)). Also, the continuum opacity sources 
and the EOS have undergone similar developments. Thanks 
to these, ID atmosphere models are in many respects very 
successful in comparisons with observations and are widely 
applied in astronomy today. 

Another approach, almost exclusively used for solar atmo- 
sphere modeUng, is the use of semi-empirical models. In these 
models, the temperature stratification is inferred from observa- 
tions (e.g. from lines forming at different heights or continuum 
center-to-limb variations). Often-used semi-empirical ID solar 
atmosphere models are the Holweger & Mueller (1974), VAL3C 
(Vemazza et al. 1976), Maltby et al. (1986) and MISS (AUende 
Prieto et al. 2001) models. A similar approach can be used 
to integrate spatially resolved observations and thus infer the 
three-dimensional (3D) atmosphere structures using inversion 
techniques (Ruiz Cobo & del Toro Iniesta 1992; Socas-Navarro 
2011). Semi-empirical modeUng is rarely attempted for other 
stars, although exceptions exist (e.g. AUende Prieto et al. 2000). 



Constructing more realistic models requires one to go be- 
yond the ID framework and model convection without relying 
on MLT. Stellar convection is an inherently 3D, time-dependent, 
non-local, and turbulent phenomenon. Therefore, one cannot 
expect ID models to reproduce all observed properties accu- 
rately, even with access to free parameters to tweak. The next 
natural step is to abandon some of these crude simplifications by 
constructing realistic 3D atmosphere models of solar convection. 
Early hydrodynamic simulations (Nordlund 1982; Nordlund & 
Dravins 1990; Steffen et al. 1989) revealed that stellar surface 
convection operates in a distinctly different fashion from the 
MLT picture. Instead of the homogeneous convective elements, 
they displayed highly asymmetrical motions with slow broad 
steady upflows interspersed with fast narrow turbulent down- 
drafts, sometimes even supersonic (e.g. Stein & Nordlund 1998, 
hereafter SN98; Asplund et al. 2000; Nordlund et al. 2009; 
Carlsson et al. 2004; Ludwig et al. 1999). The advent of 3D 
simulations, which are constructed from first principles, has 
enabled astronomers to predict various observables such as solar 
granulation properties and spectral line profiles astonishingly 
well. More recent solar 3D simulations are remarkably good at 
reproducing the observed center-to-limb variation (e.g. Pereira 
et al. 2009a; Asplund et al. 2009; Ludwig 2006). 

3D atmosphere models are by design free from the ad- 
justable parameters of MLT and other parameters such as micro- 
and macro-turbulence that have hampered stellar spectroscopy 
for many decades. Instead, in 3D simulations, convection 
emerges naturally, by solving the time-dependent hydrodynamic 
equations for mass-, momentum- and energy-conservation, cou- 
pled with the 3D radiative transfer equation in order to account 
correctly for the interaction between the radiation field and 
the plasma. Also, the non-thermal macroscopic velocity fields 
associated with convective motions are rendered realistically, 
and various natural kinetic consequences such as overshooting 
and excitation of waves emerge from the simulations, without 
the need for further ad hoc modeling or additional free param- 
eters. The inhomogeneities in the convective motions arise 
spontaneously and self-organize naturally to form a distinct 
flow pattern that exhibits the characteristic granulation at the 
surface. Furthermore, additional spectral observables such as 
limb-darkening and detailed spectral line shapes, including 
asymmetries and shifts, are also modeled unprecedentedly 
accurately with 3D models for the Sun (Nordlund et al. 2009; 
Pereira et al. 2009b) 

For metal-poor late-type stars it has been shown (Asplund 
et al. 1999b; Collet et al. 2006, 2007) that flie assumption of pure 
radiative equilibrium in the convectively stable photospheric 
layers of classical hydrostatic models is generally insufficient. 
In particular in the upper photosphere, the thermal balance is 
instead primarily regulated by radiative heating due to spectral 
line re-absorption of the continuum-radiation from below and 
adiabatic cooling due to the expansion of upflowing gas. In 
metal-poor stars, the balance between heating by radiation and 
cooUng by mechanical expansion of the gas occurs at lower 
temperatures because of the weakness and scarcity of spectral 
lines at low metallicities. By contrast, ID MLT models have no 
velocity fields outside their convection zones, and are therefore 
in pure radiative equilibrium. The temperature stratification 
there is therefore regulated solely by radiative heating and 
cooUng, thus neglecting altogether the adiabatic cooling com- 
ponent. This results in an overestimation of the temperatures by 
up to ~ lOOOK in ID models at very low metallicities, which 
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can potentially lead to severe systematic errors in abundance 
determiiiations based on ID models (see Asplund et al. 1999b; 
Asplund & Garcia Perez 2001; Ludwig et al. 2010; Collet et al. 
2009; Gonzalez Hernandez et al. 2010). These shortcomings 
of ID models are manifested as inconsistencies in the analysis 
of observed spectra, such as abundance trends with excitation 
potential of the lines (e.g. analysis of NH lines in the very 
metal-poor star HE1327-2326, by Frebel et al. (2008) and 
discrepant abundances between atomic and molecular lines 
involving the same elements (e.g. Nissen et al. 2002). For 
further discussion, we refer to a review of possible impacts 
of 3D models on stellar abundance analysis by Asplund et al. 
(2005). 

Additionally, there are discrepancies between observations 
and predictions from ID models of the solar structure in the 
context of helioseismology, which point to mistakes in the outer 
layers of theoretical ID stellar-structure models, and which are 
usually referred to as surface effects (Rosenthal et al. 1999). 
With classical ID stellar structures, higher frequency p-modes 
of the Sun are systematically shifted due to discrepancies at 
the upper turning points of the modes, which occur in the 
superadiabatic peak at the top of the convection envelope. 
Rosenthal et al. (1999) found better agreement of stellar 
structures with heUoseismic observations, when including the 
mean stratification of solar 3D models at the top, since the 
turbulent pressure, usually neglected in ID models, extends the 
resonant cavity. Also, it was found that with 3D solar models 
the predicted p-mode excitation rates are much closer to helio- 
seismic observations(Nordlund & Stein 2001; Stein & Nordlund 
2001). Ludwig et al. (2009b) compared the power spectra of 
the photometric micro-variability induced by granulation and 
found good agreement between the theoretical predictions of 
3D models and CoRoT observations (Mathur et al. 201 1). 

With the aid of 3D simulations, stellar radii have been 
derived for a number of red giants from interferometric obser- 
vations Chiavassa et al. (2010, 2012). The determined stellar 
radii are slightly larger than estimated with the use of ID 
models, which has an impact on the zero point of the effective 
temperature scale derived by interferometry. Furthermore, 
Chiavassa et al. (2012) showed that a detailed knowledge of 
the granulation pattern of planet-hosting stars is crucial for the 
detection and characterization of exoplanets. 

Several 3D magnetohydrodynamics codes with realistic 
treatment of radiative transfer have been developed and applied 
to the modeling of stellar surface convection. Here, we make 
use of the SxAGOER-code, which is developed specifically to 
run efficiently on the massively parallel machines available 
today (Nordlund & Galsgaard 1995^; Kritsuk et al. 2011). 
The BiFRosT-code is an Oslo derivative of the SxAGGER-code 
(Gudiksen et al. 2011), tailored for simulations of the solar 
photosphere and chromosphere, and therefore including true 
scattering (Hayek et al. 2010). Other widely used codes are 
CO^BOLD (Freytag et al. 2012), MURaM (Vogler et al. 
2005) and ANTARES (Muthsam et al. 2010), which have been 
independently developed in the last decades. Beeck et al. (2012) 
compared solar models from three of the above 3D stellar 
atmosphere codes (Stagger, CO^BOLD and MURaM), and 
showed that the models are overall very similar, despite the 
distinct numerical approaches. Most of the available 3D stellar 
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convection codes are now highly parallelized, which when 
coupled with the computational power available today makes it 
feasible to construct grids of 3D convection simulations within a 
reasonable time-scale. Grids of 2D and 3D atmosphere models 
already exist (Ludwig et al. 1999, 2009a; Trampedach 2007). 
Clearly, the age of 3D atmosphere modeling has arrived, partly 
driven by the rising demand created by improved high-resolution 
spectroscopic and asteroseismic observations. 

In this paper, we present a new large grid of 3D model atmo- 
spheres for late-type stars, covering an extensive range of efi'ec- 
tive temperatures, surface gravities, and metallicities. In Section 
2 we describe the methods that we followed in order to com- 
pute the 3D model atmospheres, with emphasis on the convec- 
tion code we used (Sect. 2.1) and on the tools we developed 
for scaling the grid models and post-process the results of the 
numerical calculations (Sect. 2.3). In Section 3, we present an 
overview of general properties (Sect. 3.1) of our simulations, 
and discuss the temporally and spatially averaged atmospheres 
(Sect. 3.2) from the 3D model atmospheres. We also compare 
our results with theoretical ID models in Sect. 3.3 correspond- 
ing to the same stellar parameters^. These have been computed 
with a specifically newly developed ID code that employs ex- 
actly the same EOS and opacities as the 3D simulations. Also, 
we compare our 3D model atmospheres with ID models from 
grids widely adopted by the astronomical community (MARCS 
and ATLAS). Finally, in Section 4, we summarize our findings 
and outline a roadmap for our future ambitions on the many pos- 
sible applications of the SxAOGER-grid. 

2. Methods 

2.1. The SiAGGER-code 

The 3D model atmospheres presented here were constructed 
with a custom version of the SxAGGER-code, a state-of-the-art, 
multipurpose, radiative-magnetohydrodynamics (R-MHD) code 
originally developed by Nordlund & Galsgaard (1995), and con- 
tinuously improved over the years by its user community. In 
pure radiation-hydrodynamics mode, the SxAGOER-code solves 
the time-dependent hydrodynamic equations for the conserva- 
tion of mass (Eq. 1), momentum (Eq. 2), and energy (Eq. 3) in 
a compressible flow 

dtp = -v-(pv), (1) 

d,p\ = -V ■ip\\ + T)-Vp+pg, (2) 
d,e = -V-(ev)-/?V-v + ^rad + ?visc, (3) 
coupled to the radiation field via the heating and cooling term 

fed = 47:p j'^K^iJx-S^)dA, (4) 

which is computed from the solution of the radiative transfer 
equation 

n-Vh = gK^(S;i-h) (5) 

to account properly for the energy exchange between matter and 
radiation (here p denotes the density, v the velocity field, p the 
thermodynamic pressure, e the intemal energy per unit volume^, 

^ We refer always to stellar atmospheric parameters. 

^ In the following, we will indicate the intemal energy per unit mass 

with e = e/p . 
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g the gravity, T the viscous stress tensor, ^yisc = llijTijdVi/drj the 
viscous dissipation rate, the monochromatic opacity, the 
monochromatic intensity, 7a = 1 /47r J^/^iiQ the monochromatic 
mean intensity averaged over the entrre soHd angle, and S a the 
source function). We have ignored magnetic fields in the present 
grid of 3D convection simulations, however, we will study their 
efi'ects in a future work. 

2.1 .1 . Details on the numerics 

The SxAGGER-code uses a sixth-order explicit finite-difference 
scheme for numerical derivatives and the corresponding fifth- 
order interpolation scheme. The solution of the hydrodynamic 
equations is advanced in time using an expUcit third-order 
Runge-Kutta integration method (Williamson 1980). The code 
operates on a staggered, Eulerian, rectangular mesh: the thermo- 
dynamic variables, density and internal energy per volume, are 
cell-centered, while momentum components are defined at cell 
faces. Also, in the MHD mode, the components of the magnetic 
field B (electric field E) are defined at the cell faces (edges). 
This configuration allows for a flux-conservative formulation of 
the magnetohydrodynamic equations, at the same time ensur- 
ing that the magnetic field remains divergence-free. The solu- 
tion of the discretized equations is stabilized by hyper- viscosity 
which aims at minimizing the impact of numerical diffusion on 
the simulated flow, while providing the necessary diffusion for 
large-eddy simulations with finite-difference schemes (see also 
Nordlund & Galsgaard 1995 for further details). The values of 
the numerical viscosity parameters^ are empirically tuned for the 
solar surface-convection simulation: they are set large enough to 
stabilize the numerical solution of the hydrodynamic equations 
and, at the same time, kept small enough to reduce their smooth- 
ing of the flow's structures. The same optimized values of the 
parameters are then applied to all other simulations in the grid. 

The version of the SxAGGER-code we used for this work is 
fully MPI-parallel. The parallehzation scales well with the num- 
ber of cores. For this project, the simulations were typically run 
on 64 cores. 

2.1.2. Geometrical properties 

The setup of the simulations is of the so-called box-in-a-star 
type: the domain of the simulations is limited to a small rep- 
resentative volume located around the stellar photosphere and 
including the top portion of the stellar convective envelope. The 
boundary conditions of the simulation box are periodic in the 
horizontal directions and open vertically. Gravity is assumed to 
be constant over the whole extent of the box, neglecting spheric- 
ity effects. However, since the size of the simulation domains 
correspond to only a fraction of the total radii of the stars (0.4 % 
and ~ 10% of the stellar radius for the solar simulation and for 
a typical logg -1.5 red giant simulation, respectively) such ef- 
fects can be regarded as small for the purposes of the current grid 
of models. Also, for simphcity, the effects of stellar rotation and 
associated Coriolis forces are neglected in the present simulation 
setup, as it would add two more dimensions to the grid. 

At the bottom, the inflowing material has a constant value of 
specific entropy per unit mass, which ultimately determines the 
emerging effective temperature. While the domains of our sim- 
ulations cover only a small fraction of the convective zone, the 
box-in-a-star setup is still vaUd because the bulk up-flows at the 

^ The actual values we used are ni = 0.005 and = 0.8 (see Eq. 9 in 
Kritsuketal. 2011). 



bottom boundary of the simulations carry essentially the same 
entropy value as in deeper layers and are mostly unaffected by 
entrainment with cooler downflows. At the beginning of each 
simulation, the entropy of the inflowing gas at the bottom is ad- 
justed in order to yield the desired T^g and, after that, is kept 
unchanged during the entire run (see Sect. 2.3). Furthermore, 
pressure is assumed to be constant over the whole bottom layer. 

The physical dimensions in the horizontal directions are 
chosen to be large enough to cover an area corresponding to 
about ten granular cells. The vertical dimensions are extended 
enough for the simulations to cover at least the range of -5.0 < 
logTRoss < +6.0 in terms of Rosseland optical depth (in fact they 
range on average from -7.3 < logTRoss < +7.5), which typi- 
cally corresponds to approximately six orders of magnitude in 
pressure (about 14 pressure scale heights). All of the simula- 
tions have a mesh resolution of 240^, since a resolution of about 
200^-250^ was found to be adequate (see Asplund et al. 2000). 
Five layers at the bottom and the top are reserved for the so- 
called ghost-zones: these extra layers serve to enforce boundary 
conditions for the high-order derivatives in the vertical direction. 
The spacing between cells in the horizontal direction {Ax, Ay) is 
constant, ranging from about 6 km in dwarfs to about 25 Mm in 
giants, while it varies smoothly with depth in the vertical direc- 
tion, in order to resolve the steep temperature gradients near the 
optical surface. These are the layers from where the continuum 
radiation escapes; they are characterized by a sharp transition be- 
tween stellar interior and outer layers in terms of thermodynamic 
quantities such as temperature, internal energy, and entropy that 
marks the beginning of the photosphere. Also, the steepest tem- 
perature gradients are found in the superadiabatic region j ust be- 
low the optical surface (0.0 < logTRoss < 2.0). Therefore, it is 
very important that the thin transition layer around the optical 
surface is well-resolved in order to ensure an accurate modeling 
of the radiative transfer and to avoid spurious numerical artifacts 
from insufficient spatial resolution. 



2.1.3. Equation of state 



We use the realistic equation of state (EOS) by Mihalas et al. 
(1988), which explicitly treats excitation to all bound states of 
all ionization stages, of all included elements. We have cus- 
tom computed tables for a mix of the 17 most abundant el- 
ements (H,He,C,N,0,Ne,Na,Mg,Al,Si,S,Ar,K,Ca,Cr,Fe and 
Ni). The only molecules that are included in the EOS are H2 
and , and they are treated on equal footing with the atoms and 
ions. For the solar abundances we employed the latest chem- 
ical composition by Asplund et al. (2009), which is based on 
a solar simulation performed with the same code and atomic 
physics as presented here. Our choice for the EOS, is supported 
by Di Mauro et al. (2002) who showed that solar models based 
on the EOS by Mihalas et al. (1988) show better agreement with 
hehoseismology in the outer 20Mm (> 0.97/?©), compared to 
models based on the OPAL-EOS. We inverted the Mihalas et al. 
(1988) EOS tables, hence the temperatures and the thermody- 
namic pressures are tabulated as a function of density and inter- 
nal energy. This inversion exploits the analytical derivatives pro- 
vided in the EOS tables to minimize losses in accuracy. These 
analytical derivatives are also used in the bi-cubic spUne inter- 
polation in the inverted tables. 
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2.1.4. Opacity 

We use the continuum absorption and scattering coefficients 
listed in detail and with references by Hayek et al. (2010). These 
include the sophisticated calculations by Nahar (2004)^ for the 
first three ions of all metals we include, except for K and Cr. 
These calculations are improvements over those forming the ba- 
sis for the OP opacities^ (Badnell et al. 2005). The line opacity 
is suppUed by the opacity sampUng (OS) data that was also used 
for the newest MARCS grid of stellar atmospheres Gustafsson 
et al. (2008). These are in turn based on the VALD-2 database^ 
(Stempels et al. 2001) of atomic fines and the SCAN-base^ (J0r- 
gensen 1997) of molecular lines. 

2.1.5. Radiative transfer 

The radiative heating and cooling rate (Eq. 4) is evaluated by 
solving the radiative transfer equation 



3 — - Ja-'^a, 



(6) 



where tx = ^ pKxds denotes the monochromatic optical depth 
along a given direction s, with a method similar to that by 
Feautrier (1964). The equation is solved at each time step on 
long characteristics crossing all grid-points at the simulation sur- 
face, along the vertical direction and along eight additional in- 
clined angles (two yu = cos0 and four ^-angles) by tilting the 
3D cube. Given the opacity kx and the source function S a, the 
monochromatic intensity can be obtained by solving Eq. (6) 
and the radiative heating and cooling rate computed by integrat- 
ing pKxih- S a) over sofid angle and wavelength. We use the 
Radau quadrature to determine the optimal ray directions to ap- 
proximate the angular integral in the calculation of the radiative 
heating and cooling rate as a weighted sum. For the radiative 
transfer calculations, we employ opacities as described above 
(Sect. 2.1.4). 

Computing the full monochromatic solution to the radiative 
transfer equation in 3D at each time step is extremely expen- 
sive. The cost of the radiative transfer calculations however can 
be reduced enormously by opting instead for an approximated 
solution based on the opacity binning or multi-group method 
(Nordlund 1982; Skartfien 2000). Following this method, we 
sort all sampled wavelength points into different bins based on 
the spectral range they belong to and on their associated opacity 
strength or, better, their formation depth, i.e. the Rosseland opti- 
cal depth tross (ta =1). where the monochromatic optical depth 
equals unity. In this way, wavelength points characterized by 
similar formation heights and belonging to the same spectral in- 
terval are grouped together (see Fig. 3). For each simulation, we 
use the ID temporal and spatial mean stratification to estimate 
the formation heights of the various wavelengths and sort the 
wavelengths into the different opacity bins. The bin selection 
and wavelength sorting process is performed twice during the 
simulation's relaxation phase after updating the individual mean 
stratifications, but is kept unchanged during the production runs 
that make up the time-series presented in this work. 

To each bin, we assign a mean opacity Ki which accounts 
for the contribution from both continuum and line opacities. To 
compute the mean opacities, we differentiate between diffusion 



^ http : / /www . astronomy . ohio- state . edu/ ~nahar/. 

* http : //cdsweb . u- strasbg . f r/topbase/. 

^ http://www.astro.uu. se/~vald/php/vald/. 

^ http : //www . astro . ku . dk/~uf £eg j/. 



and free-streaming limits, i.e. between the optical thick and op- 
tical thin regimes, below and above the photospheric transition 
zone, respectively. The mean bin-opacity Ki is calculated as a 
Rosseland-fike average 



J Mi) dl / Jx(i) 



dBx 
Kx dT 



dA 



(7) 



in the optical thick regime, and as a mean-intensity-weighted 
mean opacity 



Kj,i= I KxJxdA I JxdA 

JAii) I J Ail) 



(8) 



in the optical thin regime, where A(i) is the set of wavelength 
points assigned to bin /. For bin the transition from one regime 
to the other around that bins optical surface is achieved by means 
of an exponential bridging of the two averages: 



-2troi 



•Kr 



-2troi 



(9) 



All simulations presented here have been run with the ra- 
diative transfer in the strict LTE approximation, i.e. under 
the assumption that he monochromatic source function 5 a (in 
Eq. 4) is the Planck function at the local gas temperature, i.e. 
S x(T) = Bx (T). For each bin ;, we compute an integrated source 
function by summing up the contributions from all wavelength 
points in the bin; 



Si = B, 



•-I ■ 

Jao) 



BxdA 



(10) 



Collet et al. (2011) showed that, with this opacity binning 
implementation, the approximation of strict LTE results in a tem- 
perature stratification very similar to the case, where scattering 
is properly treated, as long as the contribution of scattering to 
the extinction is excluded in the calculation of the streaming- 
regime mean opacities /cj_, (Eq. 8). They also showed that in- 
cluding scattering as true absorption leads to erroneous atmo- 
sphere structures due to overestimated radiative heating in the 
optically thin layers. However, these findings have so far being 
verified only a small sample of stellar parameters, therefore we 
cannot rule out that scattering needs to be accounted for properly 
in certain cases. Nonetheless, evaluating the radiative transfer in 
strict LTE greatly eases the computational burden compared to 
the full scattering case (Hayek et al. 2010). 

The radiative transfer equation is solved for the individual 
opacity bins (Eq. 5) for all layers that have max (tross) < 300, 
while in the deeper layers, we use instead the diffusion approxi- 
mation, which is fulfilled to a high degree at such depths. With 
the opacity binning approximation, the radiative heating rate 
term (Eq. 4) takes then the form 



Anp^KiUi- 



Bi) 



(11) 



where 7, is the mean intensity computed from the solution of the 
radiative transfer equation for bin i. 

For the relaxation phase of the simulation runs we considered 
six bins, while for the final models we used twelve opacity bins. 
We have developed an algorithm for the bin selection, which 
wiU be explained further below (see Sect. 2.3.2). Towards lower 
surface gravities {\ogg < 2.0) and higher effective temperatures, 
numerical artifacts in the radiative transfer can occasionally de- 
velop and manifest as a Moire pattem in the integrated outgoing 
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Fig. 1. Kiel diagram (rejj -logg diagram) siiowing tlie targeted 
SxAGGER-grid parameters for the 217 models, comprising seven differ- 
ent metallicities (colored circles). Four additional standard stars (see 
text) are also indicated (squares). In the background, the evolutionary 
tracks for stellar masses from 0.7 to 1.5 M© and for solar metallicity are 
shown (thin grey lines). 

intensities due to very steep temperature gradients in the photo- 
sphere. For those situations, we solve the radiative transfer equa- 
tion on an adaptive mesh with finer vertical resolution, which 
is dynamically optimized to resolve regions where temperature 
gradients are steeper. The radiative heating and cooling rates 
computed on the adaptive mesh are then interpolated back to the 
coarser hydrodynamic depth scale under the consideration of en- 
ergy conservation. 

2.2. The SiAGGER-grid models 

The SxAGGER-grid covers a broad range in stellar parameters 
with 217 models in total. The range in effective temperature is 
from Teff = 4000 K to 7000 K in steps of 500 K, while the gravity 
ranges from \ogg = 1.5 to 5.0 in steps of 0.5. The grid also 
covers a broad range in metallicity starting from [Fe/H] - -4.0 
to -1-0.5 in steps of 1.0 below -1.0, and steps of 0.5 above 
that^. We decided to apply the same parameters Tgg and log^ 
for all metallicities, in order to facilitate the interpolation of 
(averaged) models within a regular grid in stellar parameters. 
In addition, the grid also includes the Sun with its non-solar 
metallicity analogs, and four additional standard stars, namely 
HD 84937, HD 140283, HD 122563 and G 64-12 that are 
presented in Bergemann et al. (2012). For metal-poor chemical 
compositions with [Fe/H] < -1.0 we applied an or-enhancement 
of [a/Fe] = -1-0.4 dex, in order to account for the enrichment by 
core-collapse supemovae (Ruchti et al. 2010). 

In Figure 1, we present an overview of our simulations in 
stellar parameter space. Therein, we also show evolutionary 

^ We use the bracket notation [X/U] = log(Nx/Nn)* - log(A'x/A'H)o 
as a measure of the relative stellar to solar abundance of element X with 
respect to hydrogen. 



tracks (Weiss & Schlattl 2008) for stars with masses from 0.7 
to 1.5M0 and solar metallicity, in order to justify our choice of 
targeted stellar parameters. Hence, the grid covers the evolution- 
ary phases from the main-sequence (MS) over the turnoff (TO) 
up to the red-giant branch (RGB) for low-mass stars. In addi- 
tion, the RGB part of the diagram in practice also covers stars 
with higher masses, since these are characterized by similar stel- 
lar atmospheric parameters. 

2.3. Scaling and relaxing 3D models 

Generating large numbers of ID atmosphere models is relatively 
cheap in terms of computational costs, but the same is not true 
for 3D models. Based on our experiences from previous simu- 
lations of individual stars, we designed a standard work-flow of 
procedures for generating our grid. More specifically, we devel- 
oped a large set of IDL-tools incorporating the various neces- 
sary steps for generating new 3D models, which we then applied 
equally to all simulations. The steps are: 

- Scale the starting model from an existing, relaxed 3D sim- 
ulation, and perform an initial run with six opacity bins, so 
that the model can adjust to the new stellar parameters. 

- Check the temporal variation of T^ff and estimate the number 
of convective cells. If necessary, adjust the horizontal sizes, 
in order to ensure that the simulation box is large enough to 
enclose at least ten granules. 

- If the optical surface has shifted upwards during the re- 
laxation, add new layers at the top of it to ensure that 

(IogTRoss>top<-6.0. 

- Determine the period ttq of the radial p-mode with the 
largest amplitude, then damp these modes with an artificial 
exponential-friction term with period ttq in the momentum 
equation (Eq. 2). 

- Let the natural oscillation mode of the simulation emerge 
again by decreasing the damping stepwise before switching 
it off completely. 

- Re-compute the opacity tables with 12 bins for the relaxed 
simulation. 

- Evolve the simulations for at least ~ 7 periods of the fun- 
damental p-mode, roughly corresponding to ~ 2 convec- 
tive turnover times, typically, a few thousand time-steps, of 
which 100 - 150 snapshots equally spaced were stored and 
used for analysis. 

During these steps the main quantities of interest are the time 
evolution of effective temperature, p-mode oscillations, and 
drifts in the values of the mean energy per unit mass and of the 
mean density at the bottom boundary, which indicate the level 
of relaxation. When the drifts in these above properties stop, we 
regard the simulation as relaxed. If these conditions were not 
fulfilled, we continued running the model, to give the simulation 
more time to properly adjust towards its new quasi-stationary 
equilibrium state. Also, when the resulting effective temperature 
of an otherwise relaxed simulation deviated more than 100 
K from the targeted T^ff, we re-scaled the simulation to the 
targeted value of T^ff and started over from the top of our list of 
relaxation steps. 

The interplay between EOS, opacities, radiative transfer and 
convection can shift the new location of the photosphere, when 
the initial guess made by our scaling procedure slightly misses it. 
This is the case for a few red giant models leading to upwards- 
shifts of the optical surface and of the entire upper atmosphere 
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Fig. 2. In the top panel, we display the non-equidistant vertical spac- 
ing Az of the depth scale as a function of geometrical depth in our solar 
model (solid line). The z-scale is optimized to resolve the flows and 
thermal structure, which naturally results in the highest spatial resolu- 
tion around the photosphere. Furthermore, we also show the maximum 
of the vertical gradient of the absorption coefficient max(dlnQ-Ross/rfz) 
as a function of depth (dotted line). In the bottom panel, we show the 
aspect ratio Az/Ax (solid line) and we also indicated its lower allowed 
limit with 1 : 4 (dashed line). The actual vertical-to-horizontal aspect 
ratio ranges from 0.26 at the photosphere to 1.18 at the bottom of the 
simulation domain. 



during the adjustment phase after the scaHng, with the average 
Rosseland optical depth ending up to be larger than required, i.e. 
{logTRoss)top ^ -6.0. In order to rectify this, we extended those 
simulations at the top by adding extra layers on the top, until the 
top layers fulfilled our requirements of {logrRoss)top < -6.0. 



The resulting new z-scale is usually not smooth, therefore 
we generate a new z-scale, which is optimized to resolve the 
region with the steepest (temperature) gradients, as shown 
in Fig. 2. The density-, energy-, and velocity cubes are 
then interpolated to this new geometrical depth scale. The 
new z-scale is constructed using the variation with depth of 
the (smoothed) maximum of the derivative of the Rosseland 
absorption coefficient, max((ilnQ'Ross/<^z), as a guide. The basic 
idea behind this approach is to vertically distribute the mesh 
points as evenly as possible on the optical-depth scale. With 
such an optimized z-scale we can efficiently resolve the same 
features with fewer grid-points, compared to an equidistant 
vertical mesh. Furthermore, a limiting vertical- to-horizontal 
aspect ratio (Az/Ax and Az/Ay) of I : 4 over the whole ver- 
tical extent is enforced. We find that this value represents 
in practice an optimal lower limit to the aspect ratio, with 
respect to numerical stability and accuracy of the solution of the 
radiative transfer equation along inclined rays. Finally, the posi- 
tion of the zero-point in the depth scale is adjusted to coincide 



with the position of the mean optical surface, i.e. (tross)(-= 



■0) '■ 



1. 



At fixed surface gravity and metallicity, the mean diameter 
of the granules, which is used for determining the horizontal ex- 
tent of the simulation, increases with higher effective tempera- 
ture (see Figs. 11 and Sect. 3.1.5). The number of granules 
present in the simulation box is retrieved with the aid of the con- 
tour routine in IDL. Based on the map of the temperature be- 
low the surface (the vertical velocity would serve equally well), 
a contour chart of the significantly hotter granules is extracted, 
from which the number of granules is counted. Concerning the 
temporal resolution of the simulation sequences of the final pro- 
duction runs, the frequency, at which snapshots are stored, is 
based on the sound-crossing time of one pressure scale height. 
Hp, in the photosphere, i.e. 



Af, 



snap - {Hp/Cs}(j=2/3) 

(see Cols. 14 and 15 in Table C.l). With the help of functional 
fits of the dependence of granule sizes and sound-crossing time 
scales on stellar parameters, the horizontal sizes of the simu- 
lation boxes and the snapshot sampling times can be estimated 
rather accurately in advance (see App. B). 



(12) 



2.3.1. Scaling the initial models 

To start a new simulation, we scale an existing one with 
parameters close to the targeted ones, preferably proceeding 
along lines of constant entropy of the inflowing gas at the 
bottom in stellar parameter space (see Fig. 6). In this way, 
we find that the relaxation process is much faster In order to 
generate an initial model for a set of targeted parameters, we 
scale temperature, density, and pressure with depth-dependent 
scaling ratios derived from two ID models, with parameters 
corresponding to the current and intended 3D model. For this, 
we used specifically computed ID envelope models (MARCS 
or our own ID models, see Sect. 3.3.1), which extend to 
logTRoss > 4.0. The reference depth-scale for all models in 
the scaling process is the Rosseland optical depth above the 
photosphere and gas pressure normalized to the gas pressure at 
the optical surface below it (logTRoss > 0.0). 

After the initial scaling, we construct the geometrical 
depth scale z for the new simulation by enforcing the same 
(quasi-)hydrostatic-equilibrium condition as in the starting 
simulation, but with the newly scaled pressure and density. 



2.3.2. Selection of the opacity bins 

As we mentioned earlier, in Sect. 2.1.5, the purpose of the 
opacity-binning approximation is to reproduce the radiative 
heating and cooling rates as accurately as possible with a small 
number of opacity-bins, in order to reduce the computational 
burden. For the assignment of wavelength points to bins, we first 
compute the opacity strengths for all of the > 10^ wavelength 
points in the opacity-sampling (OS) data from the MARCS pack- 
age (Gustafsson et al. 2008). The histogram of their distribution 
as a function of wavelength (see Fig. 3) exhibits a characteristic 
"L"-shape. Shorter wavelengths (UV) require more bins to re- 
solve the wide range in opacity strength, while the lower part of 
the L-shaped distribution at longer wavelengths (optical and IR) 
calls for a better resolution in terms of wavelength. Therefore, 
we initially make a division in wavelength at Ax, between the UV 
and the optical/IR (see boundaries of bin 1, 11 and 12 in Fig. 3) 
and comprising approximately an equal number of wavelength 
points. These two regions are then in turn subdivided evenly 
into opacity bins according to the number of /l-points. By trial 
and error, we found that a binning scheme with three bins in the 
A < Ax region and eight bins for A > Ax, one of which being a 
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Fig. 3. We show the twelve opacity bins selected for the solar sim- 
ulation by plotting the opacity strength (or, more precisely, the forma- 
tion height) against wavelength for all sampled wavelength points. The 
individual bin elements are indicated by colored symbols. For clar- 
ity, we plotted only a subset of the wavelength points considered for 
the opacity binning procedure. In the background, we included the 
smoothed histogram of the opacity strength distribution (blue contour). 
This shows how the majority of /i-points are mostly concentrated close 
to the continuum-forming layers and only a smaller fraction contributes 
to lines. 



large one and comprising the stronger lines in the optical and 
IR (bin number 10 in Fig. 3) gives a good representation of the 
monochromatic radiative heating and cooling. We iterate the bin 
selection with slight differences (e.g., one additional division in 
opacity strength for the 8 bins in the lower part of the optical and 
IR) and by small adjustments, and choose the bin selection with 
the smallest relative difference between the total heating rates 
computed with opacity binning ^bin and the full monochromatic 
solution for the average stratification of the 3D simulation, 
i.e. 



max[5^bin] = 



maxl^il 



(13) 



We found that the individual selection of some of the bins 
displays a highly non-linear response to small changes. In most 
cases an even distribution was favored by the minimization. 
Naturally, our method will typically find only a local minimum 
due to the small sample of iterations instead of a true global 
minimum. However, our method is a fast, repeatable, and 
automatic selection of the opacity bins, which minimizes the 
human effort significantly, while at the same time yielding 
very satisfactory results. Moreover, the possible deviation 
from the global minimum due to our automated bin selection 
and its resulting uncertainties are anyways smaller than the 
overall uncertainties associated with the opacity binning method. 

In Fig. 4, we compare the resulting radiative heating and 
cooling rates from the monochromatic calculation against those 
from the opacity binning solution for the mean stratification of 
our solar model. The radiative heating and cooling rates from the 
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Fig. 4. Comparison of the radiative heating and cooling resulting from 
monochromatic computations q,i (filled dots) and the opacity binning 
method q^i^ (solid line) for the solar model mean stratification. In the 
top panel we show both q^^^ vs. optical depth, while in the bottom 
panel, we compare the two against each other. 



simplified opacity binning appear rather similar to those from the 
monochromatic solution, thereby supporting our approach. For 
the solar model, our algorithm finds a bin selection that is just 
slightly less accurate (max[5<7bm] = 2.78%) than an optimized 
manual bin selection (1.86 %). Incidentally, with six bins, we get 
max[5g'bin] - 3.54%. We obtain an average max[i5^bin] for all 

while with six bins we 

; = 3.0'; 



the grid models of max[5^bin] = 2.38 ' 

get max^^bin = 3.0%. We find that max[(y^bin] increases slightly 
with Jeff and [Fe/H]. 



3. Results 

The spatially and temporally averaged mean 3D stratifications 
(hereafter (3D)) from all of our 3D models will be available on- 
line. The methods we applied to average our models are explic- 
itly described in a separate paper We provide the models in our 
own format, but also in various commonly used formats suited 
for standard ID spectrum synthesis codes such as MOOG (Sne- 
den 1973), SYNTHE (Kurucz 1993) and Turbospectrum (Plez 
2008; de Laverny et al. 2012), together with routines to interpo- 
late the (3D) models to arbitrary stellar parameters. In this paper 
the discussion will be confined to global properties and mean 
stratifications only. More extensive discussions and presenta- 
tions of the wealth of details present in the data of our 3D RHD 
models will be performed systematically in subsequent papers. 
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Fig. 5. Overview of the constant entropy value of the adiabatic convec- 
tion zone, which is given by the fixed entropy of the inflowing plasma 
at the bottom, ibot> (circles) as well as the atmospheric entropy mini- 
mum, .Sinin, (squares) for two metallicities ([Fe/H] = -0.5 and 0.0, blue 
and black respectively) against Teff . The jump in entropy, As, is in- 
dicated with vertical dotted lines. Simulations with the same gravity 
are connected with functional fits for ibot and .Smi,, (solid and dashed 
lines respectively; see App. B), while similar simulations with different 
[Fe/H] are connected with short solid black lines. 



3.1. Global properties 

In Table C. 1, we have listed the stellar parameters together with 
the thermodynamic values fixed for the inflows at the bottom, i.e. 
the internal energy Sbot^ density pbot and entropy .Sbot. as well as 
important global properties for our 3D simulations. Before we 
consider the (3D) stratifications in Sect. 3.2, we briefly discuss 
some (temporally averaged) global properties. 



3.1.1. Stellar parameters 

Surface gravity and metallicity are input parameters for a simu- 
lation, while the effective temperature is a property ensuing from 
the fixed entropy of the inflowing material at the bottom ibot- We 
calculate the effective temperature from the spatially averaged 
emergent radiative energy flux Frad and the Stefan-Boltzmann 
law Teff - [Fad/crY^^, with cr being the Stefan-Boltzmann con- 
stant. In Column 1 of Table C.l we have listed the resulting 
temporally averaged Teff of our final, relaxed simulations. These 
differ somewhat from the targeted TeffS, since we do not know 
a priori, the relation between ibot and T^-g. However, the ma- 
jority of our models (72%) deviate less than 50K, and the mean 
deviation for the whole grid is ATg^ ~ 32 K. 



3.1 .2. Constant entropy of the adiabatic convection zone 

The main input parameter that has to be adjusted is ^bot, which 
has the same value as the entropy in the deep convection zone 
due to the adiabaticity of convection, i.e. ibot = -^ad- This is also 
the reason, why the results from our rather shallow boxes are 



valid. We set .Sbot by specifying a fixed value for the density and 
energy per unit mass for the inflowing material at the bottom, 
Pbot and fibot- The actual values of ebot> Pbot and ibot applied in 
our simulations are given in Table C. 1 . Furthermore, we provide 
functional fits for .Sbot (see App. B). We compute the entropy by 
integrating the first law of thermodynamics in the form 

1 / dp\ 

ds^ —Ide-pth-^j, (14) 

adding an arbitrary integration constant in order to shift the 
zero-point of the entropy to a similar value as in Ludwig et al. 
(1999). In Fig. 5, we show .s^bot against Tg^ for [Fe/H] = 0.0 and 
-0.5, as an example. The value for .s-bot increases exponentially 
with higher T^g and with lower logg, and linearly with metallic- 
ity with a moderate slope. 

In order to increase the effective temperature solely, i.e. the 
emergent radiative flux Fy^^ at the top boundary, the total energy 
flux ascending from the convection zone has to increase by that 
same amount due to conservation of energy (see Sect. 3.2.8). 
On the other hand, when we keep Teff fixed and decrease the 
surface gravity, this in turn will cause the density to decrease 
correspondingly (see Sect. 3.2.4). Therefore, to maintain the 
same energy flux, either the transported heat content (As, s) or 
the mass flux (p or v..) needs to be enhanced. When the energy 
flux is carried by a larger mass flux, then we speak of an en- 
hancement in convective efficiency^^ (see Sects. 3.2.4, 3.2.2 and 
3.2.8). When one considers ebot with stellar parameters, then it 
clearly depicts qualitatively the same characteristic changes as 
■Shot- By inserting the perfect gas law'^ in Eq. 14 one obtains 

ds k dp 

ds , 

T uniu p 



(15) 



from which one can immediately see that the entropy increases 
with internal energy, s oc e, and also increases with lower density 
s oc — Inp. 

On the other hand, an increase in metallicity leads to a higher 
entropy of the adiabat and also a larger atmospheric entropy - 
jump (see Fig. 5). Furthermore, we find increased velocities 
(and As) and decreased densities at higher [Fe/H] (see Sects. 
3.2.2 and 3.2.4), which in turn affects the convective efficiency. 
The dependence on metallicity can be unveiled with the follow- 
ing approximation. Assuming that radiative energy transport has 
taken over the convective one near the optical surface, then, in 
the radiative diffusion approximation, we can write the radiative 
flux as 



16a- dT 
3 p/CRoss dz 



eff- 



(16) 



Restricting our analysis to the depth, where T ^ T^-g at the sur- 
face, we get 



!eff ' 



1 dT 

PKRoss dz ' 



(17) 



The values for pth(p,s) and T (p,E) are given in the EOS tables 
in the covered range of log(p/|^g/cm^j) 6 [-14,1] in 57 steps and 
log(£/[erg/g]) 6 [11, 14] in 300 steps. 

' ' In the present paper, we use the term convective efficiency as quali- 
tative description as explained in the main text, which deviates from the 
MLT definition (see App. A.l). 

p = kTp/fimu, with k being the Boltzmann constant, yu the mean 
molecular weight, and m„ the atomic mass constant. 
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Fig. 6. Showing lines of constant entropy of the adiabat i'bot from 1.3 
to 4.0 X 10'^^ erg/g/K in steps of 0.3 x lO'^"* erg/g/K and for [Fe/H] = 
-3.0 and 0.0 (blue solid and black dashed lines respectively) in the Kiel 
diagram. Additionally, we show evolutional^ tracks for 0.6 to 1.5 
with solar metallicity (thin grey lines). 



while, at the same time, the temperature and density are cou- 
pled by the ideal gas law, T cc 1/p, and hydrostatic equilibrium, 
T oc -1/Hp. Now, we hold T^g fixed, i.e. Frad ~ const., and 
increase solely the metallicity, then due to the higher [Fe/H], 
the opacity /cross increases as well, since the opacity depends di- 
rectly on the metallicity. Indeed, we find this in our simulations, 
which we show in Fig. 15 (bottom panel), with a-ross taken just 
below the optical surface at tross ~ 3.0. The temperature gra- 
dient dT/dz increases (decreases) with higher [Fe/H] for cooler 
(hotter) Teff's (see top panel in Fig. 15), in order to push the ra- 
diation through (see Sect. 3.2.7). As a consequence, the density 
decreases in order to compensate for the increase in the opacity 
(see middle panel in 18), and therefore to keep Eq. 17 constant, 
as initially assumed. One could also argue that as the atmosphere 
heats up with increasing metallicity, the hotter gas expands and 
the density p decreases, and we find Aa-ross 1/Ap. The over- 
all increase in dT/dz will also lead to a larger entropy jump as 
shown in given in Fig. 18 (top panel). 

We emphasize that the dependence of both ibot and As 
with stellar parameters is quite non-trivial, since not only is it 
coupled to the changes in the total energy fluxes, but it is also 
affected by the differences in the transition from convective 
to radiative transport of energy with stellar parameters. In 
particular, the non-local radiative transfer depends non-linearly 
on the conditions present in stellar atmospheres, especially 
changes in the opacity and the EOS will strongly influence 
the radiative transfer. Additionally, .Sbot will be influenced by 
changes in the efficiency of the convective energy transport, that 
is in the convective fluxes arising from the hydrodynamics (see 
Sect. 3.2.8 for more details). 

An analytical derivation of ibot as a function of stellar pa- 
rameters is rather difficult, as explained above. Nonetheless, we 
can fit ibot with stellar parameter in a functional form based on 



the results from our simulations (see App. B). This has been 
done previously, based on 2D RHD models with solar metallic- 
ity by Ludwig et al. (1999). In Fig. 6, we show how ibot varies 
across the Kiel diagram {T^ff - \ogg diagram) for [Fe/H] = 0.0 
and -3.0. In the case of .s^ot^ our results with solar metallicity 
are in remarkable good agreement with those of Ludwig et al. 
(1999), in fact these differ only sUghtly with ^ibot - 1-2.5%, 
despite the inherent differences between 2D and 3D convection 
simulations, the adopted EOS, opacities, and radiative transfer 
treatment. We find that ibot (which depicts .Sad) lies on lines of 
constant entropy in the Kiel diagram, in fact for the solar metal- 
licity these lines with the same entropy of the adiabat run almost 
1 -diagonally. Moreover, towards lower metallicity, we find two 
significant differences for the lines of constant ibot. the first one 
being that the slopes of the lines steepen, and the second be- 
ing that the distances in the T^ff - log g plane between the lines 
decrease. The latter implies metal-poor stars feature a broader 
range in .Sbot compared to metal-rich ones. As we will see in 
Sect. 3.2, this has important consequences for the stratifications. 

3.1.3. Entropy jump 

The upflows enter the simulation box at the bottom with the con- 
stant entropy value of the adiabatic convection zone, Sbot. and as- 
cend until they reach the superadiabatic region (S AR) just below 
the visible surface, where the convective energy is converted to 
radiative energy. In the photosphere, the mean free path for the 
continuum radiation grows large enough for the gas to become 
transparent, and the overturning upflow at the surface loses its 
internal energy as photons escape and carry away entropy. Fur- 
ther above in the nearly isothermal atmosphere with an exponen- 
tially decreasing density the entropy increases again due to the 
non-LTE EOS. This leads to a minimum, imin - min [(.?)], just 
above the surface (logTRoss < 0.0) in the temporal and horizon- 
tal averaged entropy (see Fig. 24). We determined the entropy 
jump from the difference of the entropy minimum and the fixed 
entropy at the bottom, i.e. 

As = Sbot-«min- (18) 

In order to calculate s^m, we used averages of the entropy on 
constant Rosseland optical depth'^, since it is the radiation 
losses that cause the sharp changes in the thermodynamic state 
around the optical surface and, therefore, the optical depth scale 
offers a better reference frame for comparisons. The averages 
on constant geometrical height (3D), smear out and thereby 
overestimate Smin increasingly towards higher Teff due to the 
increasing level of corrugation of iso-s surfaces on the geomet- 
rical scale (see Fig. 24). The constant entropy at the bottom Sbot, 
on the other hand, is a fixed input value for each simulation. 
It is worthwhile to mention that the main contribution to the 
variation in As as a function of stellar parameters is due to .Sbot, 
since .Smin varies just slightly with stellar parameter compared to 
Sbot (see Fig. 5). 

In Fig. 5 we show Smin (dashed Unes) as well as As (dotted 
lines), and it is obvious that the minimum in entropy increases 
just slightly with increasing Teff, while the jump increases with 
the constant entropy value of the adiabatic convection zone 
quasi-exponentially at higher Teff and lower \ogg. This can 
be concluded more easily from Fig. 18 (top panel), where we 
display Sbot vs. T^ff (see also Col. 8 of Table C.l and for As and 
Smin in App. B). We note that the location of the entropy jump 

Averages on constant column mass density yield a very similar i'min- 
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Fig. 7. We compare the entropy jump As against the constant entropy 
value of the adiabatic convection zone .Sbot for all grid models (filled cir- 
cles; color-coding explained in the legend). Furthermore, we show also 
hyperbolic tangent functional fits (see Eq. B.4 and App. B) between 
^eff = 4000 and 6000K (red lines; top panel) and between log^ = 1.5 
and 4.5 (grey lines; bottom panel). 



essentially represents the boundary of stars and the jump is to 
be regarded as physically realistic, which is a result of 3D RHD 
simulations. The entropy minimum coincides with the position 
of the upper end of the SAR. A similar sharp drop occurs for 
most of the thermodynamic quantities of interest (e, T and 
Hei), whereas p and ptot display a marked change of gradient. 
Moreover, the jump in entropy is an important value, since it is 
a direct measure of the efficiency of convective energy transport 
(which is set by the four MLT parameters, especially umlt, in 
the framework of MLT, see Bohm-Vitense 1958; Henyey et al. 
1965). Towards cool dwarfs As becomes smaller, indicating 
a higher convective efficiency, while towards hotter stars the 
large entropy jumps reffect a low convective efficiency. We 
will present a calibration of (Tmlt based on the entropy jump 
in a subsequent study, as previously done by Ludwig et al. 
(1999) using 2D convection simulations. While we found .Sbot 
to be astonishingly similar to those of Ludwig et al. (1999), 
in the case of As our findings are also qualitatively similar, 
but quantitatively the differences are here much larger with 
6As ^20-70%. 

It is quite remarkable how closely the variation of As resem- 
bles the change of .Sbot with stellar parameters (see Fig. 5). Mo- 
tivated by this, we compare directly As against ibot in Fig. 7. 
We find a nice correlation between As vs. ibot- At lower ibot 
values. As seems to converge towards 0.0 (a negative jump is 
not expected, since the atmosphere is losing energy in form of 
radiation from the photosphere), while for ibot <: 1 -7, As grows 
linearly with Sbot and only a modest level of scatter. In Fig. 7, 
we color-coded the T^ff- and logg-values respectively, to show 
how the residuals depend systematically on atmospheric param- 
eters. Models with higher T^ff (bright orange dots) and higher 
logg (dark grey dots) settle along higher As and vice versa. In 



order to illustrate this better, we have fitted a set of hyperbolic 
tangent functions (see Eq. B.4), which we show also in Fig. 
7. We included functional fits between Tgff - 4000 and 6000K 
(red/orange lines in top panel) and between logg =1.5 and 4.5 
(grey lines in bottom panel). Hence, we find hotter dwarfs along 
lines at larger As, while cooler giants settle along lines at smaller 
entropy jumps. 

Interestingly, in the linear part (sbot Z 1-7) As(sbot) displays 
a rather universal slope of As/sbot ~ 0.85, while the actual 
off'set in the ordinate depends mainly on T^s and logg. Another 
interesting aspect is that Tefr shows a similar strong influence as 
logg. The latter, however, is obviously expressed in logarithmic 
scale, therefore the influence of T^g is much stronger. On the 
other hand, when one performs a similar hyperbolic tangent 
functional fit for a fixed value of [Fe/H], then As is dispersed 
around the functional fit with such a large scatter that a fit is 
rather meaningless. Therefore, in contrast to Teff and logg we 
find no systematic trends with metallicity. 

Based on the strong correlation between the entropy jump 
As and Sbot, it is of interest to investigate what other scaling re- 
lations may be manifested for other stellar properties. With As 
as an inverse measure of convective efficiency, we expect that 
in light of such scaling relations, important quantities depend- 
ing on the entropy jump will also similarly scale systematically 
with Sbot, in particular density and velocity (see Sects. 3.2.2 and 
3.2.4), and therefore also the calibrated mixing-length of a par- 
ticular MLT implementation. 

3.1.4. Emergent intensity 

While classic ID models are inherently horizontally symmetric, 
therefore lacking a visible granulation pattern, the emergent 
intensity of 3D models features inhomogeneities exhibiting rich 
details, which arise due to the presence of turbulent convective 
motions. We give an overview the emergent intensity of our 
simulations in Fig. 8. Therein we display a main-sequence 
(MS) simulation (the Sun), a tumoff (TO) simulation, a K-giant, 
and a K-dwaif model, each with four different metallicities. To 
facilitate direct comparisons among the four metallicities, we 
kept the horizontal sizes and the color scales for the continuum 
intensities fixed from [Fe/H] = 0.0 for the individual stellar 
categories (we extended the metal-poor simulations by exploit- 
ing the periodic horizontal boundary conditions). The dark 
regions depict the cold intergranular lanes, while the brighter 
areas are the hot granules. The radiation above the granules 
originate at higher geometrical heights, while for downdrafts it 
comes from much lower heights. This is because the opacity 
is highly nonlinear due to the strong temperature sensitivity of 
the H"-opacity (/^h- ~ see SN98), which is the by far the 
dominant continuum opacity source in the visible for late-type 
stellar photospheres. Since the temperature difference between 
the granules and the intergranular lanes is very large (> 10^ K), 
one observes deeper regions in the downdrafts and higher 
altitudes in the granules. 

An immediate, interesting aspect that leaps to the eye from 
the overview presented in Fig. 8 is the qualitative self-similarity 
of the granulation patterns despite the large variations in 
size-scales. The emergent intensity increases towards higher 
Teff and decreases for lower surface gravities, as expected. 
From Fig. 8, it is also clear that the granule sizes decrease 
with metallicity (due to smaller Hp, see Sect. 3.1.5; see 
also Collet et al. 2007). Also apparent is the change of in- 
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Fig. 8. We show an overview of the emergent (bolometric) intensity for a selection of stars, namely main-sequence (MS), tumoff (TO), K-giant 
K-dwarf (from left to right, respectively) at a given time instant. For each star, we show four metallicities [Fe/H] = 0.0,-1.0,-2.0 and -3.0 
(from top to bottom, respectively). To facilitate comparisons between the different metallicity of each star, the intensity scale and the horizontal 
geometrical size of the metal-poor simulations are identical to [Fe/H] = 0.0, and the individual intensity contrasts [%] are indicated in each box. 



tensity contrast with stellar parameters, as we will discuss below. 



In order to discuss the changes in the intensity, we show 
in Fig. 9 the temporally averaged histograms of the intensity 
/ normalized to their individual mean intensity 7, thereby en- 
abling a direct comparison between different stellar parameters. 
The histograms of the intensity show two components: a peak 
at lower (darker, /// < 1.0) intensities, resulting from the cool 
downdrafts, and an often broader component at higher (brighter. 



/// > 1.0) intensities, arising from the upflowing hot granules 
(these findings are qualitatively to be expected; see SN98). 

As clearly depicted in Fig. 9, the shapes of the two 
components change with stellar parameters, in particular the 
amplitudes and widths, thereby changing the overall shape. The 
two components can be clearly extracted from histograms at 
higher T^g, where the intensity contrast is increasingly enhanced 
and eventually produces a distinctly bimodal distribution, which 
is a manifestation of the hidden or naked granulation (see Fig. 
10 and Nordlund & Dravins 1990). In order to better illustrate 
this, we also included the full width at half maxima (FWHM) of 
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Fig. 9. We show the temporally averaged histograms of the (bolo- 
metric) intensity (solid lines) for various stellar parameters (the bin- 
size is 100). To ease comparisons between the different stellar parame- 
ters, we normalized the individual intensity scales with its mean value. 
Furthermore, we have indicated the respective FWHM of the intensity 
histograms, /pwHMi (dotted lines), a measure of the intensity contrast. 
Note the different ordinate scale in the top panel. The bimodal distri- 
butions are given due to the asymmetry and the large contrast in the up 
and downflows. 



the individual intensity histograms /fwhm in the figure, to show 
that it correlates with the intensity contrast. On the other hand, 
at lower T^-ff the intensity contrast decreases in general, so that 
the two components overlap, leading to a single narrower higher 
peak in the histogram, thereby becoming indistinguishable from 
each other in the histogram. In addition, we find that the individ- 
ual contribution to the intensity from upflows and downflows is 
often asymmetric, meaning that the amplitudes of the two peaks 
in the bimodal distribution are unequal (see Fig. 9). In general, 
for dwarfs, we find that the relative importance of downflows 
with respect to upflows in terms of the peak contribution to the 
intensity distribution increases with increasing Teff. However, 
we also find exceptions, e.g. at lower metallicity where the 
behavior at Teff = 4500K is actually the opposite. Also, the 
balance between upflows and downflows varies with surface 
gravity. The intensity histograms for giants are in general 
broader (higher contrast) compared to dwarfs of the same Jeff 
(see /fwhm), hence exhibiting a larger intensity contrast. For 
dwarfs at lower metallicity (right bottom panel) the bimodality 
is more pronounced and the /fwhm (contrast) is broader (higher) 
towards higher T^ff, while at lower T^g the /fwhm (contrast) 
becomes narrower (lower) compared to solar metallicity (left 
bottom panel). The latter hints at an enhancement of the effect 
of hidden or naked granulation. 

To illustrate the latter in more detail, we show in Fig. 10 the 
rms of the bolometric intensity fluctuations for [Fe/H] - 0.0 and 
-2.0, which is commonly referred as the intensity contrast 



A/„ 



1/2 



(19) 



Fig. 10. Overview of the (bolometric) intensity contrast A/i-ms against 
Tetf for [Fe/H] = -2.0 and 0.0 (blue and black respectively). Models 
with the same gravity, but different Teff are connected. The enhanced 
naked or hidden granulation at lower metallicity can be extracted in the 
larger range of intensity contrasts (see text for details). 



with 7 being the mean intensity and the number of data points 
(see Roudier & Muller 1986). It is essentially defined as the 
relative standard deviation, hence it reflects the width of the 
intensity distribution (see Fig. 9). This often measured value is 
very suitable for quantifying the range of brightness fluctuations 
due to granulation. The intensity contrast increases with higher 
Teff and lower \ogg. For our solar simulation we get an intensity 
contrast of 15 %, which is close to the one found by SN98 with 
16 % (see Col. 10 in Table C.l). 

Towards higher Tefr, we find that .Sbot, As, and the vertical 
velocity increase, as shown in Sects. 3.1.2 and 3.2.2. For 
increasingly hotter stars, the top of the convective zone, Ztop.cz, 
penetrates higher and higher above the optical surface due 
to larger vertical velocities (see Fig. 17). Additionally, at 
higher Teff (higher Sbot ™d As), the overall temperatures and 
their fluctuations also increase, implying that one observes 
increasingly higher layers, since the dominant H~-opacity, 
hence the optical depth, depends sensitively on the temperature. 
Therefore, the granulation pattern is enhanced at higher Teff, 
while on the contrary for lower T^g the granulation becomes 
less visible, since Ztop.cz recedes below the optical surface in 
the latter case (see overview in Fig. 8). This phenomenon has 
been already described by Nordlund & Dravins (1990) as naked 
granulation. 

Interestingly, in our simulations, we find that at lower 
metallicity the effect of naked and hidden granulation is more 
pronounced, in the sense that the range of contrast from cool, 
low-contrast dwarfs, to hotter, high-contrast dwarfs, is 6 1 % 
larger (from 10.9 to 17.6) for our [Fe/H] = -2.0 simulations 
than for solar metallicity (see Fig. 10). 
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Fig. 11. Overview of the granule diameter dgran derived from the 
maximum of the mean 2D spatial power spectrum of the bolometric 
intensity against Tg^g for [Fe/H] = -2.0 and 0.0 (blue and black respec- 
tively). Models with the same gravity are connected by their respective 
functional fits (solid lines; see App. B). 



The variations in the intensity and in the intensity contrast 
with stellar parameters have important ramifications for ob ser- 
vations. At the one hand, the enhanced naked or hidden granu- 
lation at lower metallicity affects the formation of spectral lines 
and the limb darkening. On the other hand, it should also lead 
to distinct signatures in the granulation background of astero- 
seismological observations and spectro-interferometric imaging. 
Here, we limit ourselves to the discussion of the global proper- 
ties of the emergent intensity patterns, and a detailed analysis 
will be performed in subsequent papers. 

3.1.5. Granule size 

The physical dimensions of the simulations boxes is^,Sy and 
in Cols. 11 and 12 of Table C.l) are selected based on the 
mean diameter of granules (see Sect. 2.3.1 and Table C.l) 
and a target of about 10 granules in the box. Additionally, we 
measured the granule sizes by calculating the 2D spatial power 
spectrum of the bolometric intensity for the time series, and 
determining its maximum from the smoothed time average (see 
SN98). This method is quite robust despite the large variations 
in gravity. In Fig. 11, we present the measured granule sizes 
ligian (given in Col. 13 in Table C.l; see also App. B), showing 
that they become larger with smaller surface gravity. Also, the 
granules of the simulations with fixed \ogg, the lowest T^ff 
are typically ~ 50% smaller compared to the simulations with 
the hottest Teff, while for the models with the lowest metallic- 
ity they are typically ~ 30 % smaller than for the metal-rich ones. 

We find a remarkable validation for the approximation of 
the maximal horizontal extent of a granule based on mass- 
conservation considerations made by SN98 (see also Nordlund 
et al. 2009). Hereafter, we denote the following relation as the 
Nordlund scaling relation (NSR). The ascending buoyant plasma 
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Fig. 12. In the top panel, we compare the granule size approximated 
with the pressure scale height df^Q (see Eq. 20) against the mean gran- 
ule diameter rfgian (same as Fig. 1 1) for all models. The individual stel- 
lar parameters are indicated (Teff, logg and [Fe/H] with red, gray, and 
blue respectively). We indicated the position of df^ic = rfgian the solid 
diagonal line. In the bottom panel, we show also the relative residuals. 
We indicated also the mean residual (dotted line), which amounts to 
~ 30%. Here, models with the same gravity are connected (solid grey 
lines). 



inside a cylindrical granule with radius r gives rise to a vertical 
mass flux with j, - [nr^]pv^. This mass flux has to deflect and 
overturn increasingly towards the top. Due to conservation of 
mass, the upflow has to drain off sideways through the edge of 
the granule within approximately one pressure scale height Hp, 
hence resulting in a horizontal mass flux ji, - [2nrHp\pvh- The 
pressure is a quantity that is scale-invariant with stellar param- 
eters, therefore the pressure scale height is preferred over the 
density scale height. Equating and jh we can solve for the 
(maximal) granular diameter, d - 2r: 



dNSR 



(20) 



We show in Fig. 12 a comparison of the granular diameters 
estimated with c/nic and from the maximum of 2D spatial power 
spectra t/gran, which is shown in Fig. 11. The astonishing tight 
correlation can solely be interpreted as clear evidence for the 
validity of the NSR. We find that the mean pressure scale height 
taken at the height of the maximum vertical rms-velocity below 
the optical surface (logTRoss ~ 2.0, see Fig. 17) gives the best 
match between c/nic and the granule sizes. Furthermore, we also 
confirm that the relevant scale-height is that of the total pressure 
scale height. Hp = p tot /pg, since we find a better agreement 
with the latter The granular diameters found from the peak of 
2D spatial power spectra are about 30 % larger than the estimate 
from Eq. 20, i.e., dgan ~ l-3t/Nic (see lower panel of Fig. 12). 
The variation of the velocity ratio v/Jv^ in the convection zone 
is rather small (v/,/v^~ 1 .0) as both are of the order of the sound 
speed, therefore the variation in Eq. 20 stems predominantly 
from Hp. In hydrostatic equilibrium the pressure scale height is 
inversely proportional to the surface gravity (Hp cc 1/g), which 
explains the strong correlation between the granular sizes and 
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log^. On the other hand, with increasing T^ff and [Fe/H], the 
pressure scale height increases slightly because of the increase 
in the ratio of pressure and density (Hp oc ptot/p)- The ratio 
actually increases even though both values decrease, since the 
density drops with height slightly more rapidly than the pressure. 

Finally, we want to mention our finding on the filling factor 
for upflows and downflows, /up and /^n respectively. We derived 
the filling factor from the sign of the velocity field in the unal- 
tered simulations on layers of constant geometrical height. Then 
we computed the mean filling factor in the convection zone, 
which yields on average for all simulations /up - 0.65 with a 
minute deviation of cr = 0.014. Therefore, we find that the mean 
filling factor is rather universal, and close to previous findings 
by SN98 with /up ~ 2/3 and /dn ~ 1/3. 

3.2. The mean atmosphere 

In the following, we want to discuss the properties of the 
mean stratifications and the temporal and spatial averages of 
various important quantities. Unless specified otherwise, the 
(3D) stratifications presented here are averages on surfaces 
of constant Rosseland optical depth, i.e. (3D) = (3D)ross. 
Whenever we employ alternative averages in the text, e.g., on 
constant geometrical height (3D),, we indicate that explicitly. 
We remark briefly that only the averages on constant geomet- 
rical height (3D), strictly fulfill the equations of conservation 
(Eqs. 1, 2 and 3), therefore also the hydrostatic equilibrium, 
while all other averages exhibit slight deviations. This and the 
actual methods we used to compute the mean stratifications are 
discussed in a separate paper (see Magic et al. in preparation). 
For the sake of clarity, we display here only a subsample of 
our grid models including MS and RGB stars (logg = 4.5 and 
2.0, which we refer to as dwarfs and giants, respectively) with 
solar and sub-solar metallicity ([Fe/H] = 0.0 and -2.0, solar and 
metal-poor, respectively). Whenever possible, we compare with 
corresponding ID models that are obtained with our ID code 
(see App. A). 

Before continuing our discussion, we would Uke to point out 
the importance of the superadiabatic region (SAR), as it will 
be referred repeatedly in the following. It is the region, where 
the transport of energy changes character, from convective to 
radiative. The top of the SAR, where Vjad = 0, marks the top 
the convection zone, since it is the uppermost point, where the 
Schwarzschild criterion is fulfilled. At the location of the peak 
of the superadiabatic gradient, one also finds the largest fluctu- 
ations and inhomogeneities in the thermodynamic variables due 
to the non-adiabatic transition to the photosphere. Furthermore, 
it is here in the SAR, where the entropy jump and the peak in 
the vertical velocity occur. In fact, the SAR effectively repre- 
sents the physical outer boundary of the convective envelope. It 
is the most dynamic part in the interior of late-type stars, where 
the largest fluctuations are found. This is the reason why hy- 
drostatic ID modeling has the greatest challenges in this rather 
small region. 

3.2.1. Temperature stratification 

We first consider the temperature stratifications in optical depth, 
which we show in Fig. 13. We also show the corresponding 
stratifications of ID theoretical model atmosphere based on our 
ID code (dotted lines) with identical EOS and opacity tables 
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Fig. 13. (3D) stratifications of the temperature vs. optical depth for 
various stellar parameters (solid lines) and ID models (dashed lines). 
Furthermore, we have marked the positions of entropy minimum (plus), 
vertical peak velocity (diamond), and maximum in Vjajj = V - V^d (tri- 
angle). 



as the 3D models. In the continuum forming layers around the 
optical surface (-1.0 < logTRoss < 0.5), the differences between 
(3D) models with different TeffS, but same logg and [Fe/H], are 
rather small besides the shift in the temperature stratification 
corresponding to the difiference in effective temperature AFeff, 
which is to be expected since Teff ^T(t = 2/3). Well above and 
below the optical surface, on the other hand, we find significant 
diff'erences between the (3D) models depending on the stellar 
parameters. 

In the upper layers (logTRoss < -2.0) of atmospheres with 
solar metallicity, we find that the behavior in mean temperature 
is similar between 3D and ID models. On the other hand, the 
metal-poor (3D) models exhibit significantly cooler temperature 
stratifications compared to (3D) models with solar metallicity 
(Ar/r (logTRoss = -0.5) ~ -1% and ~ -14% for [Fe/H] = -1.0 
and -3.0 respectively), in particular for dwarfs {\ogg - 4.5). 
The temperature stratification in the upper photospheres of 
solar-metallicity models is largely controlled by radiative equi- 
librium, while for low-metallicity models this is not generally 
the case: for metal-poor models, the absorption features become 
considerably weaker, therefore, the radiative heating by spectral 
line re-absorption (^lad) is dominated by the adiabatic cooling 
due to expansion of the ascending gas (-/?thV ■ v) in the energy 
balance (Eq. 3), leading to an equilibrium structure at cooler 
temperatures (Asplund et al. 1999b). For cool, metal-poor 
giants (e.g., Teff = 4000 K, \ogg = 2.0), we recognize the effects 
of molecule formation on the structure of the high atmosphere. 
At sufficiently low temperatures, molecules start to form, 
which contribute with a large line opacity, shifting the balance 
from adiabatic to radiative heating and cooling, resulting in 
a stratification closer to the radiative equilibrium one (see 
Gustafsson et al. 2008). These effects are rather non-linear. 
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Fig. 14. {3D) temperature stratifications with the variation of one 
stellar parameter at a time, while the two others are fixed (Tgfi, logg, 
and [Fe/H], from top to bottom, respectively). We show our standard 
averages on constant optical depth (3D) (solid line) and on constant 
geometrical depth {3D}^ (dotted line). 
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Fig. 15. Overview of the maximum temperature gradient Vpeak (top 
panel) and Rosseland opacity a'ros.s taken at the height tros, « 3.0 (bot- 
tom panel) against Teff for [Fe/H] = -2.0 and 0.0 (blue and black re- 
spectively). Models with the same surface gravity are connected by 
their respective functional fits in the top panel (solid lines; see App. B). 



and we find no obvious systematic trends within our grid models. 

We would now like to examine the influence of individual 
stellar parameters on the temperature stratification. Therefore, 
we show in Fig. 14 the temperature stratifications of models 
where we separately vary one at the time {T^ff, logg, and 
[Fe/H]), while keeping the other two parameters constant. Fig- 
ure 14 (top panel) shows, as expected, that with increasing T^ff, 
the temperature stratification becomes overall hotter above, but 
also below the optical surface, in order to provide the required 
total energy flux (higher enthalpy, Eq. 24). We find in our 
simulations that the increased TeffS with hotter stratifications are 
accompanied by lower densities and higher vertical velocities 
below the surface (see /^itot (log tross = 2.0) in Fig. 22 represen- 
tatively for the p). The net effect on the convective flows are 
lower mass fluxes for higher TeffS, since the decrease in density 
is predominating the increase in velocity, therefore resulting in 
a more inefficient convection. This is compensated with higher 
entropy jumps (see Fig. 18 with As as an inverse measure for 
convective efficiency), hence higher temperatures and steeper 
temperature gradients. On the other hand, the temperatures in 
the upper, radiative layers increase less with increasing T^ff than 
in the deeper, convective ones. We find with decreasing surface 
gravity (middle panel in Fig. 14) the same correlations as with 
increasing TgffS before, the temperature stratifications become 
hotter below the photosphere, and due to lower densities we 
find a more inefficient convection, while the upper atmosphere 



is less affected. For lower metallicities (bottom panel), the 
temperature stratifications are significantly cooler, both above 
and below the optical surface (Ar/r(logTRoss = 2.0) ~ -5% and 
-15% for [Fe/H] = -1.0 and -3.0 respectively). At the top the 
stratifications are cooler at lower [Fe/H] due to the dominance 
of adiabatic cooling over radiative heating. Below the optical 
surface, we find higher densities with lower velocities and 
entropy jumps (while the mass flux is increasing), therefore, 
leading to an efficient convection with shallow temperature 
gradients at lower metallicities. We find cooler models that fall 
below the opacity edge, which we describe below (compare 
Teff - 5000K in Fig. 16), follow an adiabatic temperature 
stratification even in the atmosphere, which coincides with 
the rather sudden change between [Fe/H] = -1.0 and -2.0 
in Fig. 14 (bottom panel). Besides our standard averages on 
constant Rosseland optical depth, we show also the averages 
on constant geometrical depth scale (3D), (here is z fixed and 
tross = (i"Ross)z), which are systematically different, in particular 
below the optical surface, but behave qualitatively in a similar 
way with stellar parameters. 

In the sub-photospheric region (logTRoss > 0.5), 
where convection dominates, the temperature gradients 
V - d\nT /dlnptot^'^ become increasingly steeper with higher 



V increases, if only the thermodynamic pressure is included, neglect- 
ing the turbulent component. 
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Teff, reflecting the hotter interior stratifications. This can 
be illustrated with the maximum in temperature gradient, 
Vpeak = niaxV, which we show in Fig. 15 (functional fits 
are also given in App. B). The increase of Vpoak with Tgff 
is close to linear, but it seems to saturate at higher T^g (see 
Teff > 6500 K). We find that the maximum in temperature 
gradient Vpeak reproduces qualitatively a similar behavior as the 
intensity contrast A/nns with stellar parameter (compare Figs. 
10 and 15), which is consistent with the strong temperature 
sensitivity of the H" opacity. Furthermore, our metal-poor 
simulations exhibit a larger range of Vpcak-values than their 
solar metaUicity counterparts, and Vpeak is similarly enhanced 
at lower metaUicity (compare also r(logTRoss - 2.0) for dwarfs 
in Fig. 13), as the intensity contrast (see Sect. 3.1.4). Our 
cool metal-poor simulations have flatter and hot metal-poor 
simulations have steeper temperature stratifications than the 
metal-rich part of our grid (see Fig. 14). Curiously, Vpeak is 
close to constant with metaUicity for the solar T^s and log^. 

We identify three main reasons for the given variations in 
temperature gradients with stellar parameters in the SAR, which 
are rooted in the hydrodynamics and the radiative transfer: ve- 
locity field, convective efficiency, and radiative back-warming. 

1. As we discussed above (see Sect. 3.1.3), the entropy jump 
As increases with eff"ective temperature according to a power 
law (see Fig. 18). This behavior arises due to the variations 
in the radiative losses (see Sect. 3.2.8), which is accom- 
panied by changes in internal energy and density (see Fig. 
16). The velocities rise rapidly, as exhibited by the growth 

of Vj_^s and also p^^^ with T^s (see Figs. 17 and 21 respec- 
tively). Similar to Vpeak, both v^^^ and p^^^^^ occur in the 
SAR, and both increase towards higher Teff and lower \ogg. 

2. The mass mixing length am changes with stellar parameters 
(Trampedach & Stein 2011). The latter is evaluated as the in- 
verse gradient of the vertical mass flux, separately in the up- 
or downflows, hence - d\n j^^up/d\np, with j,,up being 
the vertical mass flux in the upflows. Therefore, the mass 
mixing length is composed of the gradients of the density 
and the vertical velocity, i.e. am ^ l/dlnp + l/Jlnv^. We 
find that am increases for lower As, v^^^^ and Vpeak- We will 
pubUsh our findings on the mass mixing length in a separate 
paper. 

3. In the lower photospheric layers, where the continuum 
forms, radiation is absorbed (blocked) by spectral lines; this 
implies that less radiative flux can be transported at the wave- 
lengths corresponding to spectral lines and, conversely, that 
more flux has to be pushed through continuum windows, 
an effect commonly referred to as line-blanketing. This in 
turn leads to a steepening of the temperature gradient and to 
additional heating of the sub-surface layers, also known as 
back-warming (see Gustafsson et al. 2008; Nordlund et al. 
2009). In ID models it is straightforward to quantify the ef- 
fect of back-warming, as done for example by Gustafsson 
et al. (2008), who found it to contribute a slight increase in 
temperatures below the surface (A7'/7'(tross = 10) ^ 5% for 
solar metaUicity stars with T^b ~ 5000K and log^ = 3.0). 
In our 3D RHD atmosphere models, line-blanketing and 
back-warming effects are also naturally included through 
our opacity-biiming method. Isolating the radiative back- 
warming effect in our 3D simulations is, however, a little 
more involved than in ID and we defer the analysis of this 
mechanism to a future paper in this series. 



The three mentioned effects are nonlinearly coupled and 
compete with each other, making it difficult to disentangle 
the individual contributions. A quantitative analysis will be 
presented in a later paper. 

We would like now to examine more closely a sample of 
{3D) models in the e-p-plane, as shown in Fig. 16, in order to 
better illustrate the variations with stellar parameters. One can 
clearly distinguish three different regimes: the adiabatic convec- 
tion zone, the photospheric transition, and the almost isothermal 
upper atmosphere. 

At the bottom boundary, sufficiently deep in the convection 
zone where entropy fluctuations become small, {6s} = 0.3%, 
the models follow closely the associated adiabats with s = ^bot 
(green fines). They deviate increasingly from their adiabats, as 
the top of the convection zone is approached. This is due to the 
entropy deficient downdrafts (cooled in the photosphere) becom- 
ing less diluted by the entropic upflows, as the optical surface is 
approached. For the ID models (blue dashed fines), the value 
of the entropy at the bottom of the stratifications is evidently 
overestimated, particularly at higher T^b, but this is because we 
haven't calibrated the omlt parameter here, and we have used a 
value of 1.5 for all models. The transition of energy transport 
from fully convective to fully radiative is clearly visible, since, 
at the optical surface, one finds a sharp isochoric (Ap ~ 0.0) drop 
in intemal energy (this is basically the enthalpy-jump Ah in Eq. 
24). The £-jump coincides with the sub-photospheric region 
(0.0 < logTRoss < 2.0), where the atmosphere starts to become 
transparent. The transition zone ends eventually at the optical 
surface (logTRoss - 0.0, marked with big squares). Above the 
optical surface (logTRoss < 0.0) the atmosphere is almost isother- 
mal (compare with the orange isotherms in Fig. 16), with expo- 
nentially decreasing density and almost constant intemal energy 
(Ae ~ 0.0). 

The entropy Sbot at the bottom grows exponentially with in- 
creasing Jeff and decreasing logg. We showed above that the 
entropy jump increases in a similar way (see Fig. 7). Here we 
find a similar behavior for the jump in internal energy ( Ae, hence 
Ah) in the photosphere. Moreover, we show in Fig. 16 the posi- 
tions of V^*^ and v^^^ located in the e-jump, and, again, we find 

both v^^s and V^^^ to scale exponentially with T^g. For v^^s 
we have indicated the amplitudes as well, which also increases 
exponentially with T^ff. All of the aforementioned positions are 
distributed rather regularly in the e-p-plane, while they are less 
so on the logTRoss-scale (see Fig. 13). The position of Smin is 
close to the optical surface and shows fittle variation in optical 
depth. 

At lower energies and densities in Fig. 16 (loge ~ 0.3 and 
logp ~ 3.0 to 0.0) we notice the effect of Hi and He: opacity 
in form of an edge in the opacity contours (log/CRoss — 5.0 to 
-1.0), since the bound-free absorption increases (more excited 
states) towards higher energy below the ionization energy, and 
they fade away again above it. Models that fall below this edge 
exhibit a rather different stratification. In particular, towards 
cool metal-poor dwarfs, i.e. lower T^g, higher logg, and lower 
[Fe/H], the models more closely follow adiabats than isotherms, 
in the atmosphere. This effect of the competition between ra- 
diative and dynamic heating (see beginning of this Sect.) above 
the convection zone becomes particularly evident at lower metal- 
licity (for [Fe/H] < -2.0). However, for the ID models (blue 
dashed fines), this is obviously not the case, since these always 
follow isotherms due to the enforcement of radiative equilib- 
rium. Furthermore, the cool metal-poor (3D) models also dis- 
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Fig. 16. We show the mean internal energy against mean density for dwarf models (logg = 4.5) with [Fe/H] = 0.0 and -2.0 (left and right 
respectively). The specific isocontours for the entropy i- (green), Rosseland opacity per volume, p/CRoss. (blue) and temperature T (orange) are 

underlayed. Moreover, the positions of entropy minimum .?niin (plus), optical surface (large square), vertical peak velocity ^m^, (diamonds) and 

maximum in Vsad (triangle) are marked respectively. The amplitude of v?'J^s is indicated with horizontal bars with markings in 1 km/s. We also 
included the ID models (blue dashed lines). The range in optical depth is shown from logTRos^ = -5.0 to +5.0 for each dex (small squares). 
However, we note that our simulations boxes are much deeper ((logTRoss) ~ +7.5). 



play higher densities at the optical surface, thereby spanning 
a larger p-range for different TeffS. The stratifications of sim- 
ulations of hotter dwarfs, on the other hand, depend little on 
metallicity. For the simulations, we have only plotted the range 
log 

■!"Ross 6 [—5.0,5.0], and the top of this is reached at much 
higher densities for the metal-poor dwarfs than for the solar 
metallicity dwarfs. Therefore, the density ranges covered above 
the optical surface by the individual atmospheres is small for 

metal-poor models (min[logp] 1.0 and -2.0, respectively; 

see Fig. 16). 

3.2.2. Velocity field 

Next, we consider the velocity field in our simulations, which 
arise self-consistently from the solution of the hydrodynamic 
equations. In Fig. 17 we show the rms-velocity of the vertical 
component Vj,rms, being the flux carrying component of the con- 
vective flows, and being the broadening component of spectral 
lines at disk center. 

The buoyant uprising plasma will experience increasingly 
a decrease in density towards the photosphere, hence a strong 
density gradient, and due to mass conservation, the convective 
motions will eventually overturn. Therefore, v, rms peaks in the 
S AR around log trqss ~ 1 .5 for dwarfs and ~ 2.3 for giants. Fur- 
thermore, since in the SAR, the transition region from convec- 
tive to radiative transport of energy takes place due to decrease 
in opacity and the subsequent radiative losses, here we find the 
strongest turbulent motions concomitant with the greatest fluc- 
tuations in all thermodynamical quantities (in particular entropy. 



temperature and density, see Figs. 24, 13 and 18). Further to- 
wards the interior, v^^^ns drops as the density increases. From 
the slightly sub-photospheric maximum, velocities fall off to a 
minimum above the optical surface, then increases again in the 
higher atmosphere (see Fig. 17). Towards upper layers, v'..,nns in- 
creases again due to p-modes, excited in the SAR but leaking out 
of the acoustic cavity as they have frequencies above the acous- 
tic cut-off. These traveling waves have smaller amplitudes in our 
metal poor simulations, probably due to the higher densities giv- 
ing higher mode masses (see Fig. 18). The declining velocity 
above the surface is due to the fact that the convective motions 
overshoot well above the top of the convection zone. We find the 

velocity minimum to occur between logTRoss 2.3 for dwarfs 

and -1.5 for giants. 

As expected, the magnitude of the velocity field is enhanced 
towards higher Tgff, lower logg, and higher [Fe/H], similar as 
A.?. The symmetry of the velocity profile changes with logg 
and metallicity, while it is little affected by Teff. For lower 
logg, the peak in the velocity field is increasingly shifted to 
optically deeper layers (e.g. at solar metallicity the average 
peak position for dwarfs is (logTRoss) ~ 1-5, while for giants 
it is ~ 2.0). The coolest metal-poor simulations display a flat- 
ter profile, and the position of the minimum is increasingly 
shifted towards higher layers, especially for extreme metal-poor 
dwarfs ([Fe/H] < -3.0), and the profile is therefore stretched and 
skewed. 

For comparison, we also show in Fig. 17 the convective 
velocity vmlt of our ID models determined by MLT. It is 
apparent that the general trends of increasing velocities with 
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Fig. 17. Vertical rms-velocity, v-j-ms! from the (3D> stratifications 
(solid lines) and convective velocity vmlt from the corresponding ID 
MLT models (dotted lines) as function of optical depth logTRoss for 
various stellar parameters. 



increasing Teff and [Fe/H] and decreasing logg, are common 
between the simulations and the ID MLT models, although 
much less pronounced in ID. Furthermore, vmlt drops rather 
sharply at the top of the convection zone (as given by the 
Schwarzschild criterion), as no overshooting is allowed for 
in our implementation of MLT. Several non-local variants of 
MLT exists, and they allow for overshooting, but none of them 
produce velocity profiles close to that of our simulations. We 
also want to mention the large asymmetry in velocities of the 
up- and downflows (SN98): in 3D simulations, the latter are 
much faster than the former (up to ~ 2 faster, in particular 
for cool dwarfs), contrary to what is normally assumed in ID 
descriptions of convection. 

The peak vertical rms-velocity, (see bottom panel of 

Fig. 18), is a good measure of the global magnitude of veloci- 
ties in the simulations. It also serves as a measure of the amount 
of turbulence present in the simulations. The actual values are 
also given in Col. 9 in Table C.l, together with a functional fit 
in App. B. The variation of v?^s with stellar parameters re- 
sembles that of the entropy jump As (compare top and bottom 
panel in Fig. 18), namely it also increases exponentially with 
higher Tgff and lower logg and linearly with [Fe/H]. An inter- 
esting aspect is the increase of v^^* and As with Tee, which are 
close to exponential, indicating a correlation between the two. 

The characteristic variations of v?^s correspond to the inverse 
variations of the density taken at the same heights as the peaks 
in K'l.rms (see middle panel in Fig. 18). This behavior arises due 
to conservation of mass (Eq. 1), which can be expressed as 



Vlnp = -Vlnv 



(21) 




for a stationary flow (5,p = 0). Of course, under this assumption, 
this equation is strictly speaking valid only locally, while we 
compare here averaged values. Despite that, we find that this 



7000 6500 6000 5500 5000 4500 4000 



Fig. 18. Overview of the entropy jump (top panel), maximal vertical 
rms-velocity below the surface vPfJ^j (middle panel) and the density at 
the same height Ppeak (bottom panel) vs. T^-g for [Fe/H] = -2.0 and 
0.0 (blue and black respectively). Models with the same gravity are 
connected with their respective functional fits (solid lines; see App. B). 
Note the inverse correlation between density and velocity. 



relation is, in general, qualitatively fulfilled across the whole 
depth of our simulations. Towards the optical surface, the den- 
sity decreases, which has to be compensated by faster velocities, 
in order to fulfill conservation of mass as well as sustain the 
energy flux. The velocity field profile results ultimately from 
the interplay between the vertical and horizontal acceleration 
due to buoyancy and overturning respectively. The latter in 
turn is set by the radiative losses that arises from the prevailing 
opacity conditions according to the thermodynamic state of the 
plasma (see Sect. 3.2.8). Furthermore, one can also reason 
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Fig. 19. Comparison of the entropy jump As (top panel) and the 
density at the optical depth of the maximal vertical velocity below the 
optical surface weighted by the entropy jump Ppe-ik/^^ (bottom panel) 

vs maximal vertical rms-velocity v^^s for all models. A linear fit is 
also indicated in both panels (red lines). 



that at a higher effective temperature, hence hotter temperature 
stratification, the density will be lower (ideal gas gives T cc 1/p; 
see also middle panel in Fig. 18), however, at the same time, 
more energy (enthalpy) has to be carried to the surface, which 
necessitates a faster flow (as is given in Eq. 25). The entropy 
jump, density, and velocity are coupled intimately with each 
other (the vertical mass flux is /- = pv-). Therefore, changes 
in one quantity imply corresponding variations in the values of 
the other quantities, and vice versa. The radiative energy losses 
at the photospheric transition generate the entropy fluctuations 
according to the prevailing opacity and the irradiation-duration, 
hence it sets the amplitude of the entropy jump As. On the 
other hand, the entropy deficient plasma with its density excess 
determines the buoyancy force, /b ~ Ap, and therefore the 
vertical velocities of the downdrafts. The downdrafts in turn 
will settle the upflows in order to deliver the required convective 
energy flux. The subtle details in the chain of causalities are 
non-trivial and beyond the scope of the present paper 

Similar to our finding in Sect. 3.1.3 of a scaling relation 
between the entropy jump and the constant entropy value of the 
adiabatic convection zone Ai(ibot)i we find here again another 
interesting, tight scaling relations between ppenk, v^lms ^nd As, 
which we show in Fig. 19. The values are plotted on a double 
logarithmic scale, to more clearly illustrate the power-law 
character of the relations. From the above discussion, it follows 
that the vertical velocity is also coiTelated with the constant 
entropy value ibot of the adiabatic convection zone and the 
density. We also show linear fits of the density ppeak and entropy 
jump As as a function of v?^s log-log scale (red lines in 
Fig. 19), exhibiting the slopes of log A^/log v?^s ~ 0.46 and 
logppeak/log v'?n^s ~ -1.20, hence a scaling with the respective 



slopes. 

In 3D RHD simulations, the non-thermal, macroscopic 
velocity fields arising from convective instabilities are computed 
self-consistently from first piinciples and therefore have an 
immediate physical meaning. They represent the buoyant 
motions associated with convection and its turbulent features, 
and their statistical properties carry equally important physical 
information as the mean temperature or density stratifications. 
By contrast, in ID atmosphere modeling, a free-parameter- 
dependent velocity field v'mlt is derived for the convective flux 
in MLT. Also, for radiative transfer and spectral line formation 
calculations, two ad-hoc Gaussian velocity distributions - the 
so-called micro- and macroturbulence (^turb and Xtmh, respec- 
tively) - are usually introduced to model Doppler broadening 
of spectral lines associated with non-thermal (e.g., convective, 
turbulent, oscillatory, etc.) gas motions in stellar atmospheres. 
The values of the micro- and macroturbulence parameters 
are determined by comparing synthetic and observed spectral 
line profiles and line strengths. Usually, a depth-independent 
value of the microturbulence ^turb and one global value of 
the macroturbulence Xtmb are applied in theoretical spectrum 
syntheses with ID model atmospheres. Full-3D line formation 
calculations using 3D models similar to those described here, 
have demonstrated that in late-type stars the required non- 
thermal Doppler line broadening is indeed primarily the result 
of Doppler shifts from the convective motions and to a lesser 
extent oscillations in the atmosphere (Asplund et al. 2000). As 
such this non-thermal velocity field is clearly depth-dependent, 
while micro- and macroturbulence are almost always assumed 
to be non-varying with depth. Furthermore, vmlt is solely 
assigned to satisfy the necessary amount of convective flux by 
the individual prescription of MLT. While, interestingly, vmlt 
mimics to a certain extent the run of Vj- nns in the interior for 
cooler dwarfs. We remark that this inteipretation is however 
not physically consistently motivated. Moreover, the convective 
velocity varies depending on its actual implementation (e.g. 
Bohm-Vitense 1958; Henyey et al. 1965) and, as such, vmlt 
should be interpreted and used with caution. We point out 
that one important motivation for conducting 3D RHD atmo- 
sphere models is the fact that the before mentioned spurious 
inconsistent velocities become redundant. The hydrodynamical 
simulations account consistently for only one unique velocity 
field. 

Another good measure for the turbulence of a velocity field 
is the absolute value of the vorticity |w| = |V x v|, which is shown 
in Fig. 20. The vorticity aiises below the surface in SAR due 
the overturning of the upflows and the turbulent downdrafts ex- 
periencing the density gradient. The peak in oj is associated 
with pronounced shear flows, which arise due deflection of the 
horizontal flows into downdrafts of the overturning plasma (see 
SN98). The vorticity is concentrated in tube-like structures in 
the intergranular lanes around the edges of granules. The run of 
the vorticity follows closely that of v,_i-nis (see Fig. 17). 

3.2.3. Turbulent pressure 

The turbulent pressure ptuib = pvl is a dynamical pressure term 
contributed by the plasma's vertical bulk flows. It appears when 
considering the horizontal averages of the momentum equation 
(Eq. 2), more specifically of the advection term therein. The 
ratio of turbulent to total pressure, Pturb/Ptot, shown in Fig. 21, 
follows qualitatively very closely the run of v'.._nns (compare with 
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Fig. 20. (3D> stratification of the absolute value of the vorticity vs. 
Tejf is shown for various stellar parameters. 
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Fig. 21. The fraction of turbulent pressure to total pressure Ptmb/Ptot 
vs. optical depth logTRoss is displayed for various stellar parameters. 



Fig. 17), namely, it peaks in the SAR (logrRoss ~ 1-5), reaches 

a minimum around logrRoss 2.0, and increases in the upper 

layers (a functional fit for [/?turb//'tot]peak given in App. B). 
In the SAR, the shape of the ptmb/Ptot profile with optical depth 
looks similar to a Gaussian function, however, towards lower 
Jeff and metallicity, it becomes increasingly skewed. Averages 
on constant geometrical depth (3D)- are similar, only the peak 
and the upper layers are slightly lower at higher r^ff s. 

For hotter stars, in particular metal -rich giants, the turbulent 
pressure becomes comparable to the gas pressure (pturbZ/'tot ~ 



0.4) in the SAR, and the atmosphere is increasingly supported 
by ptuib- This means that neglecting the turbulent pressure, as is 
usually done in ID models, would significantly overestimate the 
gas pressure. The consequence of this is a faulty, inconsistent 
stratification, since the overestimation in gas pressure comes at 
the cost of altering other physical quantities like the density, even 
when the temperature stratification looks similar compared to a 
(3D> model. 

We find that v,rms is very close to Vtmb = [<Pturb)/<p)]'^^, 
since the latter is basically the density-weighted analog of 
the former Therefore, one can replace the depth-independent 
^turb.iD^ which is used in ID models that include turbulent pres- 
sure (pturb.iD =ySpv2^j.j^ jjj), with v,,mis (i.e. v-,nns ~ Vturb.io)- Then 
the additional scaling factor /3 becomes redundant. Additionally, 
we want to mention that since v'tuit.iD is constant with depth, 
it will lead to an over- and underestimation of the pturb.iD above 
and below the optical surface, respectively (see Fig. 17). Clearly, 
using v,,nns instead would reduce the error We also expect v,_rms 
to correlate with the microturbulence^^ ^turb, since v,nns is a 
horizontal average of the velocity field (see also Uitenbroek & 
Criscuoli 201 1). In ID line formation calculations, ^turb is intro- 
duced in order to compensate for missing Doppler broadening in 
the line extinction profile. However, how the actual coupling of 
K-,nns with ^turb is occurs, is nontrivial due to the non-locality of 
the radiation field, which is additionally impeded by non-linear 
atomic physics. Therefore, this correlation has to be determined 
empirically by comparing the results of 3D line-formation cal- 
culations with their (3D) counterparts (see Steffen et al. 2009). 
We intend to perform such calculations in an upcoming work. 

Finally, Chiavassa et al. (2011) showed that using a realis- 
tic turbulent pressure contribution to the hydrostatic equilibrium 
in ID red supergiant atmospheres, greatly improves the derived 
surface gravity in these stars. This extra pressure component 
also leads to an expansion of the atmosphere compared to a ID 
model stratification without turbulent pressure. This is refeiTed 
to as atmospheric levitation (see Trampedach 2001). This will 
affect p-modes by affording them a larger cavity, and hence low- 
ering their frequencies. This is part of the seismic near-surface 
effect which has plagued helio- and asteroseismology. 

3.2.4. Total pressure and density 

The total pressure is defined as the sum of thermodynamic and 
turbulent pressure, ptot - Pth + Ptmb^ and the former consists of 
gas and radiation pressure, pth - Pgas + Pmd- In Fig. 22, we show 
the total pressure for various stellar parameters. In contrast to 
the previous quantities, ptot decreases with higher T^g, lower 
logg, and higher metallicity. From the three stellar parameters, 
the influence of the metallicity is the strongest. We find the 
highest pressures (and densities) in the coolest metal-poor 
dwarfs and the lowest pressures in the hottest metal-iich giants. 
In the upper layers of hot metal -poor dwarfs, we find pressures 
systematically increased with respect to their ID counterparts, 
which is accompanied by similar behavior in p, /7gas and pturb- 
As we showed above, a significant fraction of the total pressure 
is contributed by turbulent pressure in the SAR and in the 
upper layers (see Fig. 21), in particular towards higher Teff and 
lower logg. Moreover, we note that the temporal and horizontal 
(3D)^-averages from our relaxed simulations are very close to 
hydrostatic equilibrium, and the turbulent pressure contributes 

This is not necessarily the case for the macroturbulence Xtwb, which 
compensates for the missing large-scale motions that alter the shape of 
the emergent spectral line profile, but not its strength (equivalent width). 
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Fig. 22. We display the (3D> total pressure ptot against the optical 
depth logTRoss for various stellar parameters (solid lines). For compar- 
ison we also plot the ID models with dashed lines. Note the different 
ordinate scale in the top panel. 
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Fig. 23. We present the (3D> stratifications of the electron number 
density pei against optical depth logTRos, for various stellar parameters 
(solid lines). For comparison, we also show the corresponding stratifi- 
cations from ID models (dashed lines). 



significantly to this equilibrium. 

Since the mean density stratifications look qualitatively sim- 
ilar to the total pressure ones, we refrain from showing them. 
Instead we prefer to show, in Fig. 18, the peak density'^ Ppeak, 
which is the density at the height of the maximum rms-vertical 
velocity (see Fig. 17). The density ppeak increases with lower 
Tgff, higher log ^, and lower [Fe/H]. These variations with stellar 
parameters arise due to the radiative transfer, since the cooling 
and heating rates (see Eq. 4) depend on density p and opacity 
K. Furthermore, in Sect. 3.1.2, we showed the interdependence 
of density p and opacity k based on the radiative diffusion ap- 
proximation (see Eq. 16). From this relation, one can see that an 
increased opacity will lead to a lower density, kcc l/p (e.g. due 
to lower log^ or higher [Fe/H]). 

3.2.5. Electron number density 

Next, we discuss the properties of the electron number density 
Pel (Fig- 23). The electron number density drops by about 
~ 2 dex at the transition from the interior to the photosphere. This 
is due to the fact that the density itself decreases here, and due 
to the recombination of hydrogen at the photospheric transition. 
The convective flux consists to ~ 1/3 of Fion (see Sect. 3.2.8), 
therefore, as the hot ionized plasma reaches the surface, it ra- 
diates away energy, recombines, and overturns into downdrafts, 
thereby reducing the number of free electrons. The electron den- 
sity increases with higher T^ff, lower logg, and higher [Fe/H]. 
The electron pressure pei = PeiksT follows similar trends as the 
electron density in terms of variations with stellar parameters 
and depth. 



The total pressure would lead to a very similar plot. 
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Fig. 24. {3D) and (3D)j. mean stratifications (solid and dashed respec- 
tively) of the entropy i- vs. the total pressure normalized to the pressure 
at the optical surface logptot/Psurf for various stellar parameters. We 
show also the .s-stratifications of the ID models (dotted lines). Note the 
different ordinate scale in the top panel. 



3.2.6. Entropy 

Local, box-in-a-star, 3D RHD atmosphere models have well de- 
fined boundary conditions at the bottom boundary because of 
the adiabaticity of the convection zone, even though they are rel- 
atively shallow and comprise only a small fraction of the con- 
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Fig. 25. We show the (SD)^^^^ (solid lines) and (3D>. (dashed lines) 
superadiabatic gradient Vs^d vs. optical depth logr^oss for various stel- 
lar parameters. Furthermore, for comparison, we also show the cor- 
responding values from ID models (dotted lines). Note the different 
ordinate scale in the top panels. 



vection zone. Indeed, the specific entropy per unit mass of the 
plasma stays constant across most of the convective zone, in par- 
ticular for the upflows. In Fig. 24, we show the average entropy. 
Below the optical surface, the entropy converges asymptotically 
against .Sbot into deeper layers, especially the averages on con- 
stant geometrical height (3D), (dashed lines). As the hot plasma 
in the granules reaches the optical surface, it becomes transpar- 
ent, thereby a large fraction of the energy is radiated away. This 
results in a decrease in entropy, until it reaches a minimum at 
the top of the convection zone (logTRoss ~ 0.0). Further up, the 
entropy then increases again due to the decoupling of the ra- 
diation and matter above the photosphere, which results in an 
almost isothermal atmosphere. The ID models (dotted lines) 
exhibit larger entropy stratifications in the deeper convection 
zone, in particular for higher T^ff, thereby overestimating the en- 
tropy jump increasingly due to the fixed mixing-length parame- 
ter aMLT with 1.5 for all stellar parameters. 



3.2.7. Superadiabatic temperature gradient 

We limit ourselves to show only the superadiabatic gradient 
Vsad = V - Vad, since it combines the important properties of 
both the total and the adiabatic temperature gradient (V and Vad 
respectively). The peak in Vgad, which looks like a skewed Gaus- 
sian function, arises solely from the total gradient V. The supera- 
diabatic gradient peaks around log tro.s.s ~ 

1.0-2.0 (see Fig. 25), 
and becomes sub-adiabatic, i.e. Vjad < 0.0, above the optical sur- 
face at logTRo.s.s < 0.0-0.5. The entropy jump correlates directly 
with the superadiabatic gradient, since Vs^d = l/cp[ds/dlnptot] 
and one can show that ds/dz - cpjHpiy - V^d)- Hence, it is no 
surprise that they exhibit similarity in the peak amplitude and po- 
sition. In particular, the peak amplitude increases with increas- 
ing Teff and \ogg (see Ai in Fig. 18; a functional fit for V^^* is 



given in App. B). The position of V^^* on the optical depth scale 
(3D)ross (triangles), hence the position of the steepest tempera- 
ture gradient, changes slightly with stellar parameters. However, 
similar to the position of v?^s, in the e-p-plane, the distribution 

of V^^^'^ is regular (see Fig. 16), namely it shifts systematically 
towards higher e and lower p with increasing Teff. 

As it is clear in Fig. 25, one finds substantial differences 
in Vsad when comparing the two (3D) stratifications with their 
ID counterparts, namely, the ID gradients exhibit distinctively 
larger amplitudes. These differences arise partly due the miss- 
ing turbulent pressure in the ID case, but do not resolve the 
discrepancies. Furthermore, we find an asymmetrically skewed 
shape towards the optical surface in the ID gradients, which is 
shared by the geometrical averages (3D),, but is not the case for 
the averages on constant Rosseland optical depth (3D)ross. A 
main reason for the shown differences between (3D)j;)oss and ID 
comes from the averaging over layers of constant tross- The un- 
derlying VadS are rather insensitive to the deviations between the 
(3D)ross and ID stratifications, so the differences arise mainly 
due to V. Between (3D), and ID the adiabatic gradients differ in 
the sub-photospheric gradient. 



3.2.8. Transport of energy 

The individual energy fluxes are quantities worthy of further 
consideration. The energy flux is conserved only on averages 
of constant geometrical height (3D),, therefore, we show and 
discuss the latter here. The total energy flux Fxox - ^^rad + ^^conv + 
Fyisc emerges from the photosphere solely in the form of radia- 
tive energy flux. The total energy flux is supplied from the con- 
vection zone by the convective energy flux, which is the sum of 
the enthalpy flux 



■Fenth - 



P 



6 y',, with (5 y, = pv, - (pv,) 



(22) 



(6jz being the horizontal fluctuations of the vertical mass flux) 
earned in the upflows and the kinetic energy flux 



1 



-pv 



(23) 



arising from the downdrafts (see SN98 and Nordlund et al. 
2009). Since the mean kinetic energy flux Ftin is negative, 
the positive enthalpy flux Fenth is the major component of the 
convective energy flux Fcom- The enthalpy flux in turn con- 
sists of the energy fluxes due to ionization Fi^n -[^~ 



thermal heat Fth 



3 Ah 
2 P 



6 y, and acoustic (sound) waves F^, 



iPthVy) - ipth) (vz)- In Fig. 26, we show the energy fluxes F^d, 
Fenth, ^kin, ^^ion, and Fth normalized to the total emergent en- 
ergy flux crT^g (for clarity, we refrain from showing Fyisc and 
^acous, since their contribution to Ftot is very small). We vary 
one stellar parameter at a time, while the other two are fixed 
(Teff, iogg, and [Fe/H], from top to bottom in Fig. 26, respec- 
tively). Just below the optical surface (0.5 < log ptot/ /'surf < 1 -0), 
both Fkin (solid lines) and Fgnth (dashed lines) increase towards 
cool metal-poor dwarfs, i.e. lower Tetr, [Fe/H] and higher logg, 
due to higher densities and velocities. The increased reduction 
of the total flux by (< 0.0) is compensated by a simulta- 
neously higher Fgnth (> 1-0). On the other hand, in deeper lay- 
ers (log/^totZ/'surf > 1-5), both converge to similar fractions for 
all stellar parameters (-0.17 and 1.14 for and Fgnth respec- 
tively). This convergence to very similar values is rather remark- 
able, and requires a detailed careful analysis. 
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Fig. 26. Behavior of the normahzed energy fluxes F/Ftot against 
the total pressure normalized to the pressure at the optical surface 
log Plot /Psurf as a function of variations in the individual stellar param- 
eters (Tgfi, logg, and [Fe/H], from top to bottom, respectively ). In 
each panel, the various curves are shown varying one of the parameters 
while keeping the other two fixed. The individual normalized compo- 
nents of Ftot are Fgnth (dashed), F]^^ (dotted), Fion (dash-triple-dotted), 
Fth (long dashed) and Fjad (solid). Averages are evaluated at constant 
geometrical height. 



The majority of the total energy flux Ftot in the convection 
zone is carried in form of ionized hydrogen'^ with Fion - 0.67, 
while thermal heat is the second most important component 
with Fth - 0.29. The acoustic energy constitutes only a small 
fraction with Facous - 0.04. SN98 found similar fractions with 
Fkin ~ -0.10 to -0.15, Fion -2/3 and F,h ~ 1 /3 for the Sun. The 
Fion and Fth fractions, which are the major constituents of the 
enthalpy flux, undergo a significant change below the surface, as 
we show in Fig. 26 for models with different stellar parameters. 
In particular, the fraction of thermal heat Fth becomes more 
significant at the cost of Fion towards cool metal-poor dwarfs. 
The thermal flux Fth reaches a maximum (up to Fth.max - 0.5) 
just below the surface, but eventually converges close to the 
above mentioned fractions in deeper layers (long dashed lines). 

In ID MLT models, the convective flux is assumed to consist 
of the enthalpy flux only, Fconv.iD (see Appendix A.l). This 
is a result of the MLT assumption of symmetric flows which 
means the kinetic energy fluxes in the up- and downflows cancel 
exactly. As remarked by Henyey et al. (1965), the details on 
^kin, ^^ion, ^^th and Facous are not at hand due to the lack of 



The given fractions are averages of all grid models. 



a self-consistent velocity field. The energy fluxes from 3D 
RHD simulations, on the other hand, arise self-consistently 
from solving the coupled equations of radiative hydrodynamics, 
without further assumptions. 

As mentioned above, the emergent total energy flux Ftot is 
carried in the convection zone mainly by the positive enthalpy 
flux 

^enth (Eq. 22). Therefore, one can approximate the convec- 
tive energy flux with the mean jump in enthalpy Ah times the 
mean vertical mass flux of the upflows below the optical surface, 
hence 

Feonv«(A/!)<pv,). (24) 

At the transition region, the enthalpy jump Ah is primarily 
caused by the strong drop in internal energy e, hence entropy s, 
and the thermodynamic pressure work is rather small (note the 
change of ptot below the surface logTRoss > 0.0 in Fig. 22), i.e. 
Ah X! TAs, where T is the temperature at the surface. By approx- 
imating T ^ Feff, one can expect the total energy flux Ftot - o"7^efr 
to depend to first order on the mean entropy jump' ^, density, and 
vertical velocity: 

crT'^f,^{As}{p){v,}. (25) 

This approximation explains the systematic variations of As, Vy 
and p with stellar parameters, which we have observed above 
(see Figs. 5, 18 and 18, respectively). The emergent radiative 
energy flux is correlated with As, p and v^, and the respective 
composition resulting from the individual contributions varies 
with stellar parameters. 

The radiative energy flux, hence the radiative heating and 
cooling rates <7iad (Eq. 4), require a different density stratifi- 
cation for different stellar parameters due to the dependence 
of opacity on thermodynamic variables, in particular k cc l/p, 
as we showed in Sect. 3.1.2 based on the radiative diffusion 
approximation (see Eq. 17). The resulting density variations 
will induce adjustments in the vertical velocity and entropy 
jump. Furthermore, we find with increasing logg and [Fe/H] 
at a fixed T^ff, the density increases, which is compensated by 
higher As and (As oc Ap~', see Eq. 15). We would like also to 
emphasize the remarkably important (non-local) influence of the 
rather thin photospheric transition region on basically the whole 
convection zone, since the entropy deficiency of the turbulent 
downdrafts are generated mainly here. The latter sets the en- 
tropy jump and the convective driving (see Nordlund et al. 2009). 

The radiative heating and cooUng rates g'tad (Eq. 4) due to 
radiative losses enter the hydrodynamic equations as a source 
and sink term in the energy equation (Eq. 3). It is the diver- 
gence of the radiative flux ^tad = -V ■ F^ad, and a large, negative 
^lad, the cooling peak, marks the transition of energy transport 
from fully convective below the optical surface to fully radia- 
tive close to the photosphere. To better illustrate the depth de- 
pendence of ^rad and the comparison among different models, 
in Fig. 27, we show the normalized cooling and heating rates 
^rad™ ~ ~^ ■ (Frad/o"7^eff)- One Can see that the amplitude of 
^rad™ decreases with higher Feff, accompanied by an increase 
in the width of the cooling peak. The position of the maximum 



For example, Ah can be determined at the top and bottom of the 
photospheric transition region (see Fig. 16). 

Here, we prefer to use As instead of directly Ah or As due to the 
adiabaticity of convection. 
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Fig. 27. We show the normalized cooling and heating rates '^f"™ 
vs. optical depth logTRo^s for various stellar parameters. The shown 
averages are retrieved on constant geometrical height (3D)^. Note the 
different ordinate scale in the top panel. 
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Fig. 28. Comparison of the temperature stratification for the (3D) and 
the ID models by showing the difference (1D)-{3D> (colored bars) 
between logTRo.ss = -5.0 and 2.0 in the Kiel-diagram. We present four 
different metallicity ([Fe/H] = -3.0,-2.0,-1.0 and 0.0). 



absolute amplitude coincides with the position of V^^^ , since 
the cooling rate (radiative loss) is setting the entropy fluctua- 
tions, hence the superadiabatic gradient (see Sect. 3.2.7). Fur- 
thermore, this location moves into upper layers for higher T^.^ 
(from logTRoss - 1.0 up to 0.2 for T^g = 4000 to 7000K respec- 
tively). On the other hand, the width of the photospheric transi- 
tion region Aph = AlogrRoss(^rad < 0) clearly widens for hotter 
Tgff, but also, in particular, for metal-poor giants (see top right 
panel in Fig. 27). While for cool dwarfs the width is typically 
Aph a: 3.0dex, for hot metal -poor giants, it reaches Apt, ^ S.Odex 
(see, e.g., model with 5000K in right top panel of Fig. 27). 

3.3. Comparison with 1D models 
3.3.1. 1D models 

A differential comparison between ID and 3D in terms of ap- 
proaches in the modeling of stellar atmospheres is of obvious 
relevance here. Therefore, we developed a plane-parallel, hy- 
drostatic, ID atmosphere code (hereafter simply referred to as 
the ID code) that is based on a similar physical treatment as the 
MARCS code with a few simplifications (see Appendix A and 
Gustafsson et al. 2008 for more details). We employ exactly the 
same EOS and opacities as in the individual 3D models, thereby 
excluding differences due to dissimilar input physics. Also, we 
applied the (3D) models as initial stratifications for the ID mod- 
els. These mean (3D) stratifications are defined on an equidis- 
tant optical depth scale from logTRoss = -5.0 to -t-5.0 in steps 
of 0.1. The well-resolved optical depth scale reduces discretiza- 
tion errors in the ID atmosphere calculations, thereby making 
the 1D-(3D) comparison more reliable. 

In Fig. 28, we show a comparison of the ID and (3D) tem- 
perature stratifications. One can immediately extract that the up- 
per layers of the atmospheres are systematically overestimated 
in the ID models by up to ~ lOOOK, in particular for metal-poor 



stars [Fe/H] < -2.0 (for solar models the maximal difference 
is ~ 500 K). In the optically thin layers of ID models, stable 
against convection, radiative equilibrium is enforced. However, 
in the upper layers of the metal-poor (3D) models, the effect 
of the non-vanishing adiabatic cooling rate is to shift the balance 
with radiative heating to lower temperatures due to a scarcity and 
weakness of spectral lines at lower metallicities (Asplund et al. 
1999a; Collet et al. 2007). On the other hand, with the ID mod- 
els, we find systematically cooler temperatures below the photo- 
sphere logTRoss - 2.0 with up to ~ lOOOK (here there is no differ- 
ence with different metallicities). Therefore, one has to keep in 
mind that, with ID atmosphere models, and for metal-poor stars 
in particular, these severe effects on the stratifications can lead 
to large systematic errors in spectroscopic abundance determina- 
tions, up to 0.5 dex or more in logarithmic abundance, depending 
on the formation height of the individual spectral lines used in 
the analysis (e.g. Asplund et al. 1999a; Asplund & Garcia Perez 
2001; Collet et al. 2006, 2007). We wiU return to this issue using 
our new grid of 3D stellar models in subsequent investigations. 

In the ID model calculations the mixing-length parameter is 
kept constant with qtmlt - 1-5, which is the commonly applied 
value (see Gustafsson et al. 2008). However, it is well-known 
that Q-MLT varies with stellar parameters (see Ludwig et al. 1999; 
Bonaca et al. 2012, Magic et al. in preparation). Therefore, we 
caution that a single fixed value will lead to severe differences in 
atmospheric stratification. The systematic deviations beneath the 
optical surface in the temperature stratification between 1 D and 
(3D) towards cool dwarfs can be interpreted as the manifestation 
of the wrong cmlt (see Fig. 28). Furthermore, pturb is neglected 
in the ID code, which affects the stratification by reducing the 
gas pressure (see Sect. 3.2.3). 
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Fig. 29. Similar as Fig. 28, however here we compare the <3D) with 
the MARCS models for [Fe/H] = -3.0,-2.0,-1.0 and 0.0. For a better 
comparison, we applied the same temperature range as given Fig. 28. 



3.3.2. MARCS and ATLAS models 

Last, we would also like to briefly compare our (3D) stratifi- 
cations with the currently widely applied MARCS and ATLAS 
models (see Fig. 29, we show only the comparison with MARCS 
modes, since the ATLAS models look qualitatively rather simi- 
lar). We find qualitatively similar deviations as with the ID mod- 
els above. At the same time, here we also have additional differ- 
ences due to the different input physics (EOS and opacities). The 
largest differences between the (3D) and ID MARCS stratifica- 
tions of metal-poor stellar atmospheres are slightly higher, with 
~ 1300K at [Fe/H] = -3.0, while for solar metallicity the tem- 
peratures are underestimated in ID by ~ 500K at mosts. Below 
the surface, the differences amount to ~ lOOOK. The ATLAS 
models are up to ~ 850 K hotter at the top and ~ lOOOK cooler 
below the surface. In both cases, the deviations at the top in- 
crease towards lower [Fe/H]. 



4. Conclusions 

We presented here a comprehensive grid of realistic, state- 
of-the-art, three-dimensional (3D), time-dependent, radiative- 
hydrodynamic (RHD) stellar atmosphere models for late-type 
stars, covering a substantial portion of stellar parameter space, 
and provided a detailed description of the approach we followed 
for the construction of models. With the aid of our realistic 3D 
RHD simulations, we are able to access and render details of 
stellar atmospheres and subsurface convection, that are out of 
reach for ID models and also inaccessible by observations. We 
presented and discussed a number of important global physical 
properties of the simulations as well as the mean stratifications 
resulting from the vast amount of data. 

The constant entropy value of the adiabatic convection zone 
has a profound influence on several aspects and properties of 
the 3D RHD simulations. In particular, we find systematic 



correlations among the constant entropy value of the adiabatic 
convection zone, the entropy jump, and the vertical velocity, 
which we interpreted as scaling relations. In addition, we find 
that the variation in intensity contrast is enhanced at lower 
metallicity. Also, we determined that the granule size scales 
basically with the pressure scale height close to the surface, 
which can be explained in the picture of what we refer to as 
Nordlund scaling relation (NSR). 

We discussed in great detail the depth-dependent temporal 
and spatial averages of various important physical quantities. 
In particular, we determined and examined various systematic 
trends in the variations of the entropy jump, the density, and 
the vertical velocity with stellar parameters. The latter can 
be discussed by regarding the changes in the transition of 
energy transport from convective to radiative at the photosphere. 
Namely, for different stellar parameters, the coupling between 
radiation and matter through the radiative transfer necessitates 
specific physical conditions due to changes in the opacity, which 
in turn alters the density. These variations in the density on 
the other hand require adjustments in the entropy jump and the 
vertical velocity. This can be illustrated under consideration of 
the total energy flux and conservation of energy. The named 
important values are coupled with each other, and these also set 
basically the general physical framework of the stellar atmo- 
sphere. The actual particular connections of these correlations 
have to be studied carefully in more detail, thus possibly leading 
to an improved understanding of the physical mechanisms 
operating in subsurface convection, hence stellar atmospheres. 

We compared our 3D models and their mean stratifications 
to ID models employing the same input physics, thereby 
revealing important systematic differences between the two 
kinds of models due to the incomplete treatment of convection 
by the ID mixing-length theory (MLT) and the assumption of 
radiative equilibrium. The latter leads to an overestimation of 
the temperature stratification in metal-poor stars. While below 
the optical surface, we find that the temperatures are typically 
underestimated due to a fixed mixing length. Also, we find that 
MLT fails to render a realistic vertical velocity field. The often 
neglected turbulent pressure has towards giants a non-negligible 
contribution on the total pressure, thereby, indicating that the 
thermal gas pressure is also overestimated significantly. We 
also quantified the differences with widely used ID atmosphere 
models, in particular ATLAS and MARCS. For a number 
of important values we provide functional fits with stellar 
parameters, so that these can be accessed immediately. Thereby, 
one can easily scale new 3D models based on these informations. 

The present work is meant to be an introduction to a series 
of papers on the SxAGGER-grid. The material discussed here is 
in fact just a small fraction of the actual information contained 
in the complete simulation data set. On the other hand, the list 
of potential applications for 3D models is also long. Because 
of space constraints, we will discuss in more detail the different 
methods we have applied for computing temporal and spatial 
averages of our 3D RHD models in a separate paper. Therein, 
we will also discuss our routines for the interpolation of (3D) 
atmosphere models, and explore the differences between (3D) 
and 3D line formation. We have also compiled details on 
statistical properties, which we will present individually. These 
can be later on utilized for a so-called 1 .5D line formation (see 
Ayres et al. 2006). An obvious plan on the agenda is also to 
analyze carefully the detailed properties of granulation and the 
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intergranular lane for different stellar parameters across the 
HR-diagram in the future. 

As the major purpose of theoretical atmosphere models, we 
are computing full 3D synthetic spectra with OPTIM3D (see 
Chiavassa et al. 2009) for all of our models (Chiavassa et al. in 
preparation). These will be made publicly accessible and can 
be used for various applications, e.g. spectroscopic parameter 
determination. Furthermore, we will derive new limb-darkening 
coefficients based on the full 3D synthetic spectra (Magic 
et al. in preparation), which are vital for applications such 
as the characterization of transiting exoplanets (Hayek et al. 
2012). We intend to calibrate photometric colors and radial 
velocities from the spectra. For abundance determinations, we 
will construct an extensive hbrary of synthetic high-resolution 
spectral lines. These will serve for deriving 3D effects on 
line-shifts, line-asymmetries and bisectors. 

Also an important application is the calibration of free 
parameters in ID models by employing our 3D simulations, 
in particular MLT and micro- and macroturbulence. We will 
incorporate our (3D) models into stellar evolutionary models, 
which will result in improved stellar structures useful for 
asteroseismology. And with 3D line formation calculations we 
will caUbrate the microturbulence. 

Despite the enormous success and the ab-initio nature of 3D 
atmosphere modeling, as last we want to mention the weak- 
nesses. In order to keep the computation costs reasonable, the 
radiative transfer is simplified with the opacity biiming method, 
which may influence the outcome. Also, the numerical resolu- 
tions of these so called large-eddy simulations are not resolv- 
ing the microscopic viscous dissipation length scales, hence, the 
need to introduce numerical diffusion. However, these do not 
affect the main properties of the macroscopic flows and of the 
physical stratification. Also, we minimized the diffusion coeffi- 
cients once under the constraint of numerical stability, and then 
applied for all the simulations. Additionally, accounting for non- 
LTE effects is extremely expensive for 3D simulations, therefore, 
these are usually neglected. An immediate rectification will be 
enhancement of the numerical resolution, in particular for giant. 
We will conduct improvements in the radiative transfer by em- 
ploying more bins and possibly more more angles in the future. 
By inserting magnetic fields with different field strength, we will 
explore the influence of the latter on convection for different type 
stars. 
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Appendix A: The STAOQER-grid 1 D atmospheres 

The following discussion concerns solely ID atmosphere mod- 
els and MLT, therefore, similar quantities as discussed above 
may deviate (e.g. Fconv)- The numerical code that we used for 
computing ID atmospheres for the SxAGGER-grid models solves 
the coupled equations of hydrostatic equilibrium and energy flux 
conservation in ID plane-parallel geometry. The ID models use 
the same EOS and opacity package in order to allow consis- 
tent 3D-1D comparisons. The set of equations and numerical 
methods employed for their solution are similar to those of the 
MARCS code (Gustafsson et al. 2008) with a few changes and 
simplifications that will be outUned in the following. The result- 
ing model atmospheres yet maintain very good agreement with 
MARCS models (see Sect. 3.3.2). 



Appendix A. 1: Basic equations 

Assuming ID plane-parallel geometry with horizontal homo- 
geneity and dominance of hydrostatic equihbrium over all verti- 
cal flow simplifies the equation of motion (Eq. 2) to the hydro- 
static equilibrium equation 



dr, 



std 



(/'gas + Ptmb) - — = 0, 



(A.1) 



where /Cstd and Tstd are a standard opacity and corresponding op- 
tical depth (e.g. the Rosseland mean), pg^s and pturb denote gas 
pressure and turbulent pressure, p is the gas density, and g is the 
surface gravity. Radiation pressure is neglected, as in the 3D 
simulations. Turbulent pressure is estimated using the expres- 
sion 



Pturb =j8pv, 



turb' 



(A.2) 



with the scaling parameter /3 that corrects for asymmetries in the 
velocity distribution and the mean turbulent velocity vturb that is 
used as a free, independent parameter. 

The depth-integral of the energy equation (Eq. 3) reduces to 
the flux conservation equation, 



^rad + 



-a-T. 



eS ■ 



■0, 



(A.3) 



where Fj-kj is the radiative energy flux, Fconv is the convective 
energy flux, cr is the Stefan-Boltzmann constant and Jeff is the 
stellar effective temperature. Contrary to the 3D case, effective 
temperature now appears as a boundary value and is thus a free 
parameter. Owing to numerical instabilities of the formulation, 



Eq. A.3 is replaced in the higher atmosphere (tross ^ 10 ^) with 
the radiative equilibrium condition 



'I 



(A.4) 



where and S ^ are the mean intensity and the source function, 
similar to Eq. (6). In the 3D case, is explicitly calculated 
and is nonzero in general. Enforcing the condition of radiative 
equilibrium q^A = in ID leads to an atmospheric stratification 
where an exact balance of radiative heating and cooling in each 
layer is achieved, ignoring the effects of gas motion. 

The mean intensity and the radiative energy flux at each 
depth are obtained by solving the radiative transfer equation. 



p-r- =h-Sx, 

dTA 



(A.5) 



where p = cos 6* with the polar angle off the vertical, I a is 
the specific intensity at wavelength A, and ta is the vertical 
monochromatic optical depth (with = above the top of the 
atmosphere). A Planck source function 5^ = is assumed. The 
monochromatic mean intensity and radiative flux are then deliv- 
ered by the integrals 



Ja 



\lhdp 

2;rJ" hpdp. 



(A.6) 
(A.7) 



In the absence of an explicit convection treatment, convec- 
tive energy transfer is estimated using a variant of the mixing 
length recipe described in Henyey et al. (1965). The convective 
flux is given by the expression 



-auLTSAcpTpvuLT, 



(A.8) 



where p is the gas density, Cp is the specific heat capacity, T 
is the temperature, and vmlt is the convective velocity. The 
weU-known free mixing length parameter aMLT = Im/Hp sets 
the distance /„, in units of the local pressure scale height Hp 
over which energy is transported convectively. See Gustafsson 
et al. (2008) for details of the expressions used to obtain the con- 
vective velocity vmlt and the factor SA = F/ (1 + F) Vsad> which 
scales super-adiabaticity Vjad = V - of the atmospheric strat- 
ification (see also Sect. 3.2.7), by a convective efficiency fac- 
tor r = cppvMLTTeiy + T~-^)/i^o'T^) with the optical thickness 
Te = KRosslm- We adopt the same parameters y = 0.076 and v = 8 
as Gustafsson et al. (2008) for the radiative heat loss term and 
turbulent viscosity that enter the above quantities. 

Appendix A.2: Numerical methods 

The system of equations is solved using a modified Newton- 
Raphson method with an initial stratification of temperature T 
and gas pressure pg^s on a fixed Rosseland optical depth grid. 
Discretized and Unearized versions of the hydrostatic equation 
and the energy flux equation (or radiative equilibrium condition, 
respectively) provide the inhomogeneous term and the elements 
of the Jacobian matrix for the system of 2N linear equations, 
where is the number of depth layers. The radiation field is 
computed for each Newton-Raphson iteration using the integral 
method, based on a second-order discretization of the fundamen- 
tal solution of the radiative transfer equation (Eq. A.5). 
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The corrections AT and A/5gas derived from the system of 
Unear equations are multiphed by a variable factor < 1 that is 
automatically regulated by the code to aid convergence. Conver- 
gence is assumed when the (relative) residuals of the 2N equa- 
tions decrease beneath a preset threshold. Note that, contrary 
to the 3D simulations, the effective temperature is now an ad- 
justable parameter; the requirement of minimal residuals auto- 
matically leads to an atmospheric stratification with correct Teff 
through the energy flux equation. 

In order to obtain a ID model, a given {3D) stratification pro- 
vides the initial input for the Newton-Raphson iterations, along 
with the targeted effective temperature and surface gravity. The 
same EOS tables that were used for the 3D simulation provide 
gas density, specific heat capacity, and adiabatic gradient as a 
function of T and p^as- Likewise, the tables containing group 
mean opacities and the Rosseland mean opacity provide the re- 
quired microphysics for solving the radiative transfer equation, 
ensuring maximal consistency with the 3D simulations. 

Once convergence has been achieved for the ID stratifica- 
tion, the mixing length parameter qtmlt can be caUbrated to ob- 
tain a better approximation to the (3D) stratification in the con- 
vection zone beneath the stellar surface. 



Appendix B: Functional fits 

The resulting amount of data from our numerical simulations is 
enormous. A convenient way to provide important key values is 
in form of functional fits, which can be easily utilized elsewhere 
(e.g. for analytical considerations). In the present paper we 
have frequently discussed various important global properties 
that are reduced to scalars. Some of them are global scalar 
values and some are determined at a specific height from the 
{3D) stratifications, i.e. temporal and spatial averages on layers 
of constant Rosseland optical depth. We fitted these scalars 
with stellar parameters for individual suitable functions, thereby 
enforcing a smooth rendering. However, we would like to warn 
against extrapolating these fits outside their range of validity, 
i.e. outside the confines of our grid. Also, one should consider 
that possible small irregularities between the grid steps might 
be neglected, which arise due to non-linear response of the EOS 
and the opacity. On the other hand, we provide also most of the 
actual shown values in Table C.l. 

We use three different functional bases for our fits and 
we perform the least-squares minimization with an automated 
Levenberg-Marquardt method. Instead of the actual stellar pa- 
rameters, we employ the following transformed coordinates: 
X = (Teff - 5777)/1000, >• = log^ -4.44 and z = [Fe/H]. Fur- 
thermore, we find that in order to accomplish an optimal fit 
with three independent variables, fiix,y,z), simultaneously, the 
metaUicity should be best included implicitly as nested functions 
in the form of second degree polynomial (z) - Tj^^QOiz', each 
resulting in three independent coefficients a,-. The linear function 

/i ix,y,z) = Uz) + xaiz)+y^ciz) (B.l) 

is applied for the following quantities: imin (Fig. 5), logppeak 
(Fig. 18), logyPr* (Fig. 18), log<igran (Fig. 10), logAfgran and 
^pcak resulting coefficients are given in Table B.l. On the 
other hand, we considered the exponential function 

fi ix,y,z) = fi ix,y,z)+^diz)cxp[x^e iz)+y^fiz)] (B.2) 



Table B.l. The coefficients a; of the hnear function /i (Eq. B.l) for ^nun 
[10"erg/g/K], logppeak [IQ-^g/cm^], log^^^, [lOkm/s], logrfgran 
[Mm] and logAfgran [lO^s]. 



«(■ 




IgPpeak 


J peak 


Ig'^gran 


IgAjgran 


ao 


1.5440 


0.3968 


-0.4626 


0.2146 


-0.7325 


fli 


0.0387 


-0.2549 


0.0568 


0.0666 


-0.0054 


ai 


0.0046 


-0.0344 


0.0068 


0.0108 


-0.0016 


bo 


0.0621 


-0.4232 


0.1988 


0.1174 


0.0410 


bi 


-0.0189 


0.1260 


-0.0255 


0.0187 


0.0046 


bi 


0.0013 


-0.0007 


0.0050 


0.0033 


0.0000 


CO 


-0.0898 


0.6814 


-0.1845 


-1.0922 


-0.9970 


C'l 


0.0038 


-0.0282 


0.0116 


-0.0462 


-0.0038 


Cl 


-0.0004 


-0.002 1 


-0.0006 


-0.0075 


-0.0006 



Table B.2. The coefficients a,- of the functional bases h and 
(Eqs. B.2 and B.3) for ibot [10"erg/g/K], As [10"erg/g/K], 

[Pturb/Ptot]peak. Vak and V^^^. 





•5bot 


As 


peak 
P turb/tot 


^peak 


ypeak 
sad 


ao 


1.5789 


-0.0006 


0.0321 


1.0941 


0.8713 


ai 


0.0455 


0.0043 


0.0459 


-0.0089 


0.0338 


ai 


0.0111 


0.0064 


0.0111 


0.0000 


0.0076 


bo 


0.0784 


0.0017 


0.0138 


0.2498 


0.3401 


bi 


-0.0183 


0.0049 


0.0007 


-0.0532 


-0.0717 


bi 


0.0071 


0.0060 


0.0019 


-0.0050 


-0.0091 


CO 


-0.1076 


0.0028 


-0.0260 


-0.4004 


-0.4847 


C\ 


-0.0028 


-0.0029 


-0.0087 


0.1052 


0.0990 


Cl 


-0.0042 


-0.0032 


-0.0016 


0.0142 


0.0155 


do 


0.1602 


0.1979 


0.1335 


-0.0600 


-0.0622 


d\ 


0.0618 


0.0675 


-0.0257 


0.0016 


-0.0006 


di 


0.0062 


0.0059 


-0.0081 


-0.0133 


-0.0128 


eo 


1.2867 


1.1479 


0.5894 






e\ 


-0.0824 


-0.0866 


0.1141 






ei 


0.0970 


0.0788 


0.0337 






/o 


-1.2136 


-1.0996 


-0.5330 






/i 


-0.0338 


-0.0316 


-0.0864 






/2 


-0.0764 


-0.0614 


-0.0249 







for Sbot, As (Figs. 5 and 18) and [pturb/Ptotlpeak (Fig- 21). For 
Vpeak and V^a*15 we applied the following function 

h {x,y,z) = fi {x,y,z) + x^Q (z) , (B.3) 

with coefficients for fi and fj, Usted in Table B.2. Finally, we 
showed in Fig. 7 the entropy jump As as a function of Sbot> which 
we fitted with 

/4 (x) — ao + a\x + aT, tanh [04 -1- as,x\ . (B.4) 
The resulting coefficients are hsted in Table B.3. 

Appendix C: Tables 

In Table C.l we have listed important global properties with 
the stellar parameters. Due to the lack of space, we show only 
an excerpt with solar metaUicity ([Fe/H] - 0.0). The complete 
table is available at CDS http : //cds . u- strasbg . f r. 
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Table C.l. The stellar parameters: effective temperature, surface gravity (Cols. 1 and 2 in [K] and [dcx]). The main input variables: the density 
Pbot> internal energy per unit mass fibot (Cols. 3 and 4 in j^lO'^g/cm-'j ,|^10^erg/gj). We also added the temperature Tbot, thermodynamic pressure 
and entropy values ibot at the bottom (Cols. 5, 6 and 7 in [K],|^10^dyne/cm^j,[l0^erg/g/Kj). Furthermore, the jump in entropy As, the 

maximal vertical rms-velocity v^_^s ^^d intensity contrast A/ms values are given (Cols. 8, 9, 10 in [lO'erg/g/K],[km/s],[%]). Finally, we 
display the horizontal and vertical box size (Cols. 1 1 and 12 in [Mm] , [Mm]), the mean granule diameter rfgran (Col. 13 in [Mm]), the time 
step Af and total time t (Cols. 14 and 15 in [lO^s] , [lO^s]). 



TeS logg IgPbot Igebot Ig^bo 



lgp_^' 



■5bot 



As 



peak 



A/n 



Igd. 



gran 



IgAf 



Igf 



4023 
4052 
3938 
4569 
4532 
4492 
4530 
4513 
4516 
4512 
4932 
5013 
4998 
5001 
4978 
4953 
4963 
5465 
5560 
5497 
5510 
5480 
5768 
6023 
5993 
5998 
6437 
6483 
6918 



1.50 
2.00 
2.50 
2.00 
2.50 
3.00 
3.50 
4.00 
4.50 
5.00 
2.00 
2.50 
3.00 
3.50 
4.00 
4.50 
5.00 
3.00 
3.50 
4.00 
4.50 
5.00 
4.44 
3.50 
4.00 
4.50 
4.00 
4.50 
4.50 



0.717 
1.125 

1.691 
0.679 
1.357 
1.785 
2.103 
2.419 
2.721 
3.013 
0.042 
0.883 
1.534 
1.960 
2.292 
2.604 
2.885 
1.084 
1.663 
2.139 
2.486 
2.791 
2.367 
1.130 
1.865 
2.301 
1.384 
2.008 
1.545 



1.124 
1.004 

0.908 
1.187 
1.060 
0.955 
0.900 
0.858 
0.835 
0.819 
1.291 
1.202 
1.082 
0.987 
0.919 
0.877 
0.854 
1.215 
1.119 
1.010 
0.947 
0.905 
0.997 
1.266 
1.122 
1.026 
1.263 
1.134 
1.283 



4.272 
4.233 
4.239 
4.342 
4.279 
4.266 
4.269 
4.277 
4.292 
4.308 
4.535 
4.374 
4.308 
4.295 
4.293 
4.301 
4.314 
4.403 
4.345 
4.322 
4.322 
4.330 
4.336 
4.493 
4.364 
4.344 
4.495 
4.386 
4.543 



1.061 
1.368 

1.889 
1.120 
1.669 
2.029 
2.322 
2.625 
2.927 
3.226 
0.700 
1.358 
1.882 
2.243 
2.538 
2.837 
3.118 
1.589 
2.062 
2.456 
2.769 
3.060 
2.688 
1.737 
2.281 
2.644 
1.989 
2.448 
2.201 



2.300 
2.018 
1.775 
2.411 
2.039 
1.808 
1.682 
1.578 
1.500 
1.434 
2.757 
2.376 
2.024 
1.805 
1.661 
1.560 
1.485 
2.337 
2.040 
1.791 
1.649 
1.547 
1.725 
2.395 
1.991 
1.771 
2.315 
1.969 
2.292 



0.602 
0.361 

0.174 
0.723 
0.395 
0.210 
0.126 
0.069 
0.037 
0.021 
1.047 
0.706 
0.399 
0.223 
0.125 
0.068 
0.038 
0.685 
0.428 
0.226 
0.128 
0.072 
0.186 
0.751 
0.397 
0.222 
0.686 
0.386 
0.673 



5.145 
4.167 

3.210 
5.845 
4.391 
3.486 
2.903 
2.367 
1.937 
1.541 
8.331 
5.880 
4.407 
3.608 
2.896 
2.363 
1.868 
5.815 
4.598 
3.597 
2.959 
2.323 
3.293 
6.183 
4.514 
3.572 
5.818 
4.516 
5.737 



18.4 

17.1 

14.4 

18.4 

17.2 

14.5 

12.2 

9.4 

7.7 

6.3 

50.4 

18.0 

16.9 

14.8 

11.7 

8.8 

6.8 

17.7 

17.4 

15.3 

12.1 

9.0 

14.6 

17.9 

17.9 

16.1 

18.3 

18.7 

18.6 



3.820 
3.322 
2.740 
3.380 
2.845 
2.342 
1.778 
1.146 
0.602 
0.146 
3.544 
2.903 
2.362 
1.813 
1.279 
0.699 
0.176 
2.447 
1.903 
1.362 
0.845 
0.301 
0.903 
1.903 
1.415 
0.845 
1.447 
0.903 
1.041 



3.490 

2.971 

2.446 

3.069 

2.517 

1.966 

1.442 

0.895 

0.399 

-0.102 

3.127 

2.586 

2.055 

1.496 

0.952 

0.455 

-0.048 

2.125 

1.572 

1.023 

0.503 

0.001 

0.601 

1.703 

1.095 

0.552 

1.221 

0.624 

0.781 



3.121 

2.623 

2.041 

2.778 

2.243 

1.643 

1.079 

0.544 

-0.000 

-0.553 

3.544 

2.204 

1.663 

1.114 

0.580 

0.000 

-0.301 

1.748 

1.204 

0.663 

0.146 

-0.398 

0.204 

1.301 

0.716 

0.146 

0.748 

0.204 

0.342 



2.176 
1.695 

1.188 

1.740 

1.241 

0.692 

0.188 

-0.319 

-0.824 

-1.301 

1.876 

1.287 

0.765 

0.220 

-0.292 

-0.824 

-1.301 

0.819 

0.318 

-0.284 

-0.770 

-1.301 

-0.678 

0.467 

-0.155 

-0.721 

-0.081 

-0.638 

-0.638 



4.352 
3.871 
3.364 
4.041 
3.282 
2.692 
2.364 
1.681 
1.276 
0.875 
4.052 
3.463 
2.765 
2.317 
1.749 
1.217 
0.796 
2.995 
2.415 
1.892 
1.230 
0.699 
1.419 
2.564 
1.942 
1.279 
2.016 
1.403 
1.362 



Table B.3. The coefficients a,- of the hyperbohc tangent function /4 (Eq. 
B.4) for fitting As as function of ibot- 





ao 


ai 


02 


03 


a4 


4000.0 


1.2910 


-0.3559 


1.0367 


-2.6408 


1.2059 


4500.0 


5.1768 


-2.1859 


4.5280 


-1.4475 


0.6756 


5000.0 


7.0730 


-3.0946 


6.8382 


-1.2330 


0.5799 


5500.0 


7.6382 


-3.4144 


7.5981 


-1.1812 


0.5636 


6000.0 


6.8963 


-2.9796 


6.9907 


-1.1769 


0.5504 


logg 




a\ 




fl3 




1.5 


5.3693 


-2.0610 


5.3770 


-1.2576 


0.5461 


2.0 


1.1012 


-0.2599 


0.9218 


-2.8316 


1.2958 


2.5 


1.5805 


-0.5023 


1.2081 


-2.5467 


1.1888 


3.0 


5.2106 


-2.1433 


4.6691 


-1.4254 


0.6548 


3.5 


4.9821 


-2.0989 


4.2522 


-1.5136 


0.7111 


4.0 


8.0957 


-3.5548 


7.9721 


-1.1979 


0.5625 


4.5 


14.1757 


-6.3782 


17.4802 


-0.8936 


0.4180 
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